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Abstract 

The review is devoted to the elucidation of the basic problems arising in the theoretical 
investigation of systems with Bose-Einstein condensate. Understanding these challenging 
problems is necessary for the correct description of Bose-condensed systems. The prin- 
cipal problems considered in the review are as follows: (i) What is the relation between 
Bose-Einstein condensation and global gauge symmetry breaking? (ii) How to resolve the 
Hohenberg-Martin dilemma of conserving versus gapless theories? (iii) How to describe 
Bose-condensed systems in strong spatially random potentials? (iv) Whether thermo- 
dynamically anomalous fluctuations in Bose systems are admissible? (v) How to create 
nonground-state condensates? Detailed answers to these questions are given in the review. 
As examples of nonequilibrium condensates, three cases are described: coherent modes, 
turbulent superfluids, and heterophase fluids. 
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1 Principal Theoretical Problems 



In recent years, the topic of Bose-Einstein condensation has been attracting very high attention. 
There have been published the books [1,2] and a number of review articles (e.g. [3-12]). This great 
attention is mainly due to a series of beautiful experiments with trapped atoms, accomplished 
in many laboratories of different countries and promising a variety of interesting applications. 
The interpretation of experiments requires the development of theory. It is well known that 
there is nothing more practical than a good theory. Only a correct theory allows for the proper 
understanding of experiments and can suggest appropriate and realistic technical applications. 

The theory of real systems with Bose-Einstein condensate was advanced by Bogolubov [13-16] 
who considered uniform weakly nonideal low-temperature Bose gas. Extensions to nonuniform 
zero -temperature weakly interacting gas are due to Gross [17-19], Ginzburg and Pitaevskii [20], 
and Pitaevskii [21]. This approach has been the main tool for describing Bose-condensed systems, 
since the majority of initial experiments with trapped atoms had been accomplished with weakly 
interacting Bose gases at low temperatures, using the techniques of cooling and trapping [22]. 

Since London [23], it is assumed that superfluidity in 4 He is accompanied by Bose-Einstein 
condensation, although detecting the condensate fraction in helium is a rather difficult experi- 
mental task. The existence in liquid helium of Bose-Einstein condensate with zero momentum 
has not been directly proved, without model assumptions, though the majority of experiments 
are in agreement with the existence of condensate fraction of about 10% [24]. The possibility 
that superfluidity is accompanied by mid-range atomic correlations [25-27] or that it is due to 
the appearance in superfluid helium of a condensate with a finite modulus of momentum [28-31] 
have also been discussed. In his works on superfluid helium, Landau [32] has never assumed 
the condensate existence. This is why the direct observation of Bose-Einstein condensation of 
trapped atoms has become so important and intensively studied phenomenon [1-12]. 

The trapped Bose gases are dilute and can be cooled down to very low temperatures. Usually, 
they also are weakly interacting. Thus, cold trapped atomic gases have become the ideal object 
for the application of the Bogolubov theory [13-16]. 

However, by employing the Feshbach resonance techniques [33,34] it is possible to vary atomic 
interactions, making them arbitrarily strong. In addition, the properties of trapped gases at 
nonzero temperature have also to be properly described. But the Bogolubov approximation, 
designed for weakly interacting low-temperature systems, cannot be applied for Bose systems at 
finite interactions and temperature. 

Attempts to use the Hartree-Fock-Bogolubov approximation resulted in the appearance of 
an unphysical gap in the spectrum [35,36]. While there should be no gap according to the 
Hugenholtz-Pines [37] and Bogolubov [16] theorems. This gap cannot be removed without loos- 
ing the self-consistency of theory, which would lead to the distortion of conservation laws and 
thermodynamic relations [16]. The situation was carefully analyzed by Hohenberg and Martin 
[38], who showed that, as soon as the global gauge symmetry, associated with the Bose-Einstein 
condensation, is broken, any theory, in the frame of the grand canonical ensemble, becomes ei- 
ther nonconserving or acquires a gap in the spectrum. This dramatic conclusion is known as 
the Hohenberg-Martin dilemma of conserving versus gapless theories. In this review, it is shown 
how a correct self-consistent theory has to be developed, being both conserving and gapless, and 
being valid for finite temperatures and arbitrary interactions. 

In the Bogolubov approach, the global gauge symmetry U(l) is broken, which yields Bose- 
Einstein condensation. Hence, this gauge symmetry breaking is a sufficient condition for conden- 
sation. But maybe it is not necessary? Some researchers state that Bose-Einstein condensation 
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does not require any symmetry breaking. This delusion is explained in the review, where it 
is emphasized that the gauge symmetry breaking is the necessary and sufficient condition for 
Bose-Einstein condensation. 

In recent literature on Bose systems, there often happens a very unfortunate mistake, when 
one omits anomalous averages, arising because of the gauge symmetry breaking. But it is straight- 
forward to show that this omission is principally wrong from the precise mathematical point of 
view. To get an excuse for the unjustified omission of anomalous averages, one ascribes such 
an omission to Popov, terming this trick "Popov approximation". Popov, however, has never 
suggested such an incorrect trick, which can be easily inferred from his original works [39,40]. 

The general self-consistent theory, presented in the review, is based on the Bogolubov shift 
of field operators, which explicitly breaks the gauge symmetry. The theory is valid for arbitrary 
interacting Bose systems, whether equilibrium or nonequilibrium, uniform or nonuniform, in the 
presence of any external potentials, and at any temperature. External potentials of a special 
type are spatially random potentials. For treating the latter, one often uses perturbation theory 
with respect to disorder. However, it is possible to show that such perturbation theory can be 
misleading, yielding wrong results. In this review, a method is described that can be used for 
disorder potentials of any strength. 

One of the most confusing problems, widely discussed in recent literature, is the occurrence 
of thermodynamically anomalous particle fluctuations in Bose-condensed systems. In the review, 
a detailed explanation is given that such anomalous fluctuations cannot arise in any real system, 
since their presence would make the system unstable, thus, precluding its very existence. The 
appearance of such anomalous fluctuations in some theoretical calculations is caused by technical 
mistakes. 

The usual Bose-Einstein condensate corresponds to the accumulation of particles on the 
ground-state level. An important problem, considered in the review, is whether it would be 
admissible to create nonground-state condensates. A positive answer is given and it is explained 
how this could be done and what would be the features of such condensates. 

Throughout the paper, the system of units is employed, where the Planck constant h = 1 and 
the Boltzmann constant ks — 1. 

2 Criteria of Bose-Einstein Condensation 
2.1 Einstein Criterion 

Bose-Einstein condensation implies macroscopic accumulation of particles on the ground-state 
level of a Bose system. This means that, if the number of condensed particles is N and the 
total number of particles in the system is N, then Bose-Einstein condensation occurs, when N 
is proportional to N. To formulate this criterion in a more precise way, it is necessary to invoke 
the notion of the thermodynamic limit, when the number of particles N, as well as the system 
volume V, tend to infinity, with their ratio remaining finite: 



iV ->■ oo , 



V ->■ oo , 



N 

— — > const . 



(1) 



Then the Einstein criterion is formulated as the limiting property 




(2) 
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where the thermodynamic limit (1) is assumed. This is a very general criterion that, however, 
does not hint on how the condensate particle number N should be found. 



2.2 Yang Criterion 

The Yang criterion [41] is related to the notion of the off-diagonal long-range order related to the 
behavior of reduced density matrices [42]. The first-order reduced density matrix p(r, r') defines 
the limit 

lim p(r, 0) = lim — , (3) 

in which r = |r|. One says that this matrix displays the off-diagonal long-range order and 
Bose-Einstein condensation occurs, when 

limp(r,0)>0. (4) 

The Yang criterion can be useful for uniform systems, but is not suitable for confined systems, 
where the limit of p(r, 0), as r — > oo, is always zero, while condensation can happen [3,9]. 

2.3 Penrose-Onsager Criterion 

Penrose and Onsager [43] showed that the occurrence of condensation is reflected in the eigenvalue 
spectrum of the single-particle density matrix. For the latter, the eigenproblem 



/ 



p(r, r')ip k (r') dr' = n k ip k (r) (5) 



defines the eigenfunctions <Pfc( r ) an d eigenvalues n k , labelled by a quantum multi-index k. The 
largest eigenvalue 

N = supn k (6) 

k 

gives the number of condensed particles Nq. That is, condensation occurs, when 

lim > 0. (7) 

This criterion is quite general and can be used for uniform as well as for nonuniform systems. 

2.4 Order Indices 

A convenient criterion can be formulated by means of the order indices for reduced density 
matrices [44-47] . Order indices can be introduced for any operators possessing a norm and trace 
[48]. Let A be such an operator. Then the operator order index is defined [48] as 

-(-4) - ™ , (8) 
log \TrA\ 

where the logarithm can be taken to any base. Considering fi\ = [p(r, r')] as a matrix with 
respect to the spatial variables results in the order index 

f ~\ - lo g||Pill (a \ 
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Using the expressions 

||pi|| = supn fe = N , Trpx = N , 

k 

yields the order index for the density matrix 

(~ \ logiVo (m\ 

u{pl) = bglv • (10) 

This order index makes it possible to give the classification of different types of order: 

^(pi) < (no order) , 

< u)(pi) < 1 (mid — range order) , (11) 

cu(pi) = 1 (long — range order) . 

The latter corresponds to Bose-Einstein condensation, when 

lim u( Pl ) = 1 , (12) 

7V-»oo 

in agreement with the previous criteria. Generally, there can exist Bose systems with mid-range 
order [45-48]. In such systems there is no Bose-Einstein condensate but there happens a quasi- 
ordered state that can be called quasicondensate [39]. 

The order indices are useful in studying confined systems. But for confined systems, the notion 
of thermodynamic limit is to be generalized. For this purpose, one has to consider extensive 
observable quantities [49,50]. Let be such an observable quantity for a system of N particles. 
The most general form of the thermodynamic limit can be given [12,51] as the limiting condition 

iV — > oo , An — > oo , — > const . (13) 

Similar conditions with respect to the system ground-state energy imply the system thermody- 
namic stability [52]. 

2.5 Condensate Existence 

The condensation criteria show that Bose-Einstein condensation imposes the following restriction 
on the behavior of the density- matrix eigenvalues n&. Recall that, by its definition, n^ means the 
particle distribution over the quantum multi-indices k. According to Eqs. (6) and (7), one has 

oc ^ -> (N ->■ oo) . (14) 



sup fc n k N 

If condensation occurs into the state labelled by the multi- index fc , so that 

sup n k = n ko , 
k 

then the condensation condition [12] is valid: 

lim — = (N — ?> oo) . (15) 

k— >ko n k 

Writing N — > oo, implies, as usual, the thermodynamic limit in one of the forms, either as in Eq. 
(1) or as in Eq. (13). 
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3 Gauge Symmetry Breaking 



3.1 Gauge Symmetry 

The global gauge symmetry U(l) for a Hamiltonian H[ijj], which is a functional of the field 
operator ip, means that this Hamiltonian is invariant under the gauge transformation 

(r) -> ^{v)e ia , (16) 

where a is a real number. That is, 

H[^e ia \ = H[1>\ . (17) 

Here and in what follows, the time dependence of field operators is assumed but is not shown 
explicitly, when it is not important and cannot lead to confusion. 
The field operator can always be decomposed into an expansion 

^(r) = ^a fcV 9 fc (r) (18) 

k 

over an orthonormal complete basis. Though, in general, the basis can be arbitrary, for what 
follows, it is important to choose the natural basis, composed of natural orbitals [42]. By defini- 
tion, the basis is natural if and only if it is composed of the eigenf unctions of the single-particle 
density matrix, defined by the eigenproblem (5). Then the eigenvalues n k describe the particle 
distribution over the quantum indices k. 

Bose-Einstein condensation can occur not to any state but only into one of the states of the 
natural basis, that is, into one of the natural orbitals. Denoting the related natural orbital by 
y?o( r ), one can write 

V>(r)=Vo(r)+Vi(r), (19) 
separating the part corresponding to condensate, 

^o(r) = cwo(r) , (20) 

from the part related to uncondensed particles, 

^i(r) = ^2a k ip k (r) . (21) 

fe^O 

By construction, the condensate part is orthogonal to that of uncondensed particles: 

J VS(r)Vi(r)dr = 0, (22) 

which follows from the orthogonality of natural orbitals. And by the definition of the natural 
orbitals as eigenfunctions of the single-particle density matrix, the quantum-number conservation 
condition 

( a l a p) = h P {a\a k ) (23) 
is valid. Because of the latter, one has the particular form of the quantum conservation condition 

(VtoVi(r)> = 0. (24) 



8 



The number-operator for condensed particles is 



N = J ^(r)^o(r) dr = aja . (25) 
And the number-operator for uncondensed particles is 

N, = ! Vl(r)Vi(r) dr = J^ a W ■ (26) 

So that the total number-operator reads as 

A> = N + A>! . (27) 
The number of condensed particles is the statistical average 

N ee (A> ) = (4a ) . (28) 
According to the condensation criteria, Bose-Einstein condensate appears when 

lim <^>0. (29) 

TV^oo N 

Till now, no symmetry breaking has been involved in the consideration. Because of this, 
one could naively think that no symmetry breaking is necessary for treating Bose condensation. 
However, the above consideration is yet nothing but a set of definitions. To understand whether 
gauge symmetry breaking is compulsory for treating Bose condensation, one has to analyze the 
properties of the defined quantities. 

3.2 Symmetry Breaking 

There are several ways how the Hamiltonian symmetry could be broken. The oldest method is by 
incorporating in the description of the system an order parameter with a prescribed properties 
corresponding to a thermodynamic phase with the broken symmetry, as is done in mean-field 
approximations [32]. Another traditional way, advanced by Bogolubov [15,16], is by adding to 
the Hamiltonian symmetry-breaking terms, getting 

^] = ^]+ e rM, (30) 

where (T[ip]) £ ex N and e is a small number. The statistical averages, with Hamiltonian (30), are 
denoted as (• • -) e . Upon calculating such an average, one should take, first, the thermodynamic 
limit N — > oo, after which the limit e — > 0. The so defined averages are called quasiaverages. It is 
also possible to combine these two limits in one, prescribing to e a dependence on iV and taking 
the sole thermodynamic limit. The latter procedure defines thermodynamic quasiaverages [53]. 
Other methods of symmetry breaking are described in the review [54]. Here, for concreteness, 
the standard way of symmetry breaking by means of the Bogolubov quasiaverages will be used. 
Spontaneous breaking of gauge symmetry occurs when 

lim lim 1 f \(Mr))e\ 2 dr > 0. (31) 
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This can also be rewritten as 

lim Km Ml^ >0 . (32) 
By the Cauchy-Schwarz inequality, 

|(a ) £ | 2 < (ala ) £ (33) 

for any e. This means that gauge symmetry breaking yields Bose condensation. 

Theorem 1. When gauge symmetry is spontaneously broken, then there exists Bose-Einstein 
condensate. 

Proof. Spontaneous breaking of gauge symmetry corresponds to Eq. (32). In view of the 
Schwarz inequality (33), it follows that 



Wo) 



lim lim x u "" /e > , (34) 



which implies Bose-Einstein condensation. 
3.3 Ginibre Theorem 

The Hamiltonian of a Bose system is a functional of the field operator ip that can always be repre- 
sented as the sum (19) of two terms (20) and (21). Thus, Hamiltonian (30) is H £ [ip] = H £ [if) ,ifji\. 
For an equilibrium system, this Hamiltonian defines the grand thermodynamic potential 

n £ = -TlnTr exp{-/3H £ [^ , ^]} , (35) 

where T is temperature and (3 = 1/T. Let us replace the operator term ip by a nonoperator 
quantity r], getting 

Q ve ee -TlnTr ex^-^fo,^]} , (36) 
and assuming that this thermodynamic potential is minimized with respect to rj, so that 

Q V£ = inf Q X£ . (37) 

X 

Ginibre [55] proved the following proposition. 

Theorem 2. In thermodynamic limit, the thermodynamic potentials (35) and (36) coincide: 

lim % = lim (38) 

This theorem holds true irrespectively from whether there is Bose condensation or not. But if, 
when minimizing potential (36), one gets a nonzero 77, then, according to condition (31), there is 
spontaneous gauge symmetry breaking. Hence, because of Theorem 1, Bose condensation occurs. 



10 



3.4 Bogolubov Theorem 

Let 

C e ^oM = {...4--M--^o...^i) e (39) 
be a class of correlation functions being the averages, with respect to the Hamiltonian H £ [ipQ, ipi], 
of the normal products of the field operators (20) and (21). And let 

C e (77,^ 1 ) = (...7 7 *...^ 1 t ...77...Vi>^ (40) 

be a class of correlation functions being the averages, with respect to the Hamiltonian H £ [rf,ipi], 
of the normal products of the field-operator terms, where the operators ipQ have been replaced by 
a nonoperator quantity i] that minimizes the thermodynamic potential (36). Then the Bogolubov 
theorem [16] holds. 

Theorem 3. In thermodynamic limit, the corresponding correlation functions from classes 
(39) and (40) coincide: 

lim C e {Mi)= C £ (rj^ 1 ). (41) 

N->oo N— >oo 

As particular consequences from this theorem, it follows that 

lim , i im m /(^o( r )^o(r)) e dr = lim ^ \r]{r)\ 2 dr , 

e^O N-^oo iV J Af->oo iV J 

lim lim (^o(r)) £ = r?(r) . (42) 
Invoking the conservation condition (24) yields 

lim lim (^i(r)) e = , 

e^O N— >oo 

lim lim (ip(r)) £ = 77(1") . (43) 

Hence, if 77 is not zero, the spontaneous gauge symmetry breaking takes place. Respectively, Bose 
condensation occurs. This important consequence of the Bogolubov theorem can be formulated 
as the following proposition. 

Theorem 4. Spontaneous gauge symmetry breaking implies Bose-Einstein condensation: 

lim lim W = i im i im faU>£. (44) 

3.5 Roepstorff Theorem 

The above theorems show that spontaneous gauge symmetry breaking is a sufficient condition 
for Bose-Einstein condensation. The fact that the former is also the necessary condition for the 
latter was, first, proved by Roepstorff [56] and recently the proof was polished in Refs. [57,58]. 

Theorem 5. Bose-Einstein condensation implies spontaneous gauge symmetry breaking: 



(aja ) 



|2 



lim < lim lim . (45) 

In the left-hand side of inequality (45), the average is taken without explicitly breaking the 
gauge symmetry. Combining theorems 4 and 5 leads to the following conclusion: 

Conclusion. Spontaneous gauge symmetry breaking is the necessary and sufficient condition 
for Bose-Einstein condensation. 
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4 General Self-Consistent Approach 



4.1 Representative Ensembles 

A statistical ensemble is a pair {J 7 , p} of the space of microstates J 7 and a statistical operator p. 
Defining the statistical operator, it is necessary to take into account all conditions and constraints 
that uniquely describe the considered statistical system. This requirement was emphasized by 
Gibbs [59,60] and Ter Haar [61,62]. Such an ensemble is termed a representative ensemble. The 
general formulation of the representative ensembles and their properties have been given in Refs. 
[54,63,64]. 

Constraints, imposed on the system, can be represented as the statistical averages of condition 
operators Ci, with i = 1,2,... being the index enumerating the condition operators. This gives 
the set of statistical conditions 

(Ci) = a . (46) 
Taking into account the latter defines the grand Hamiltonian 

H = H + Y J \C l1 (47) 

% 

in which H is the energy operator and \ are Lagrange multipliers guaranteeing the validity of 
conditions (46). 



4.2 Bogolubov Shift 

The most convenient way of gauge symmetry breaking for Bose systems is by means of the Bo- 
golubov shift [16] of the field operator, when the field operator ip of a system without condensate 
is replaced by the field operator 

^(r)=77(r)+^(r), (48) 

in which i](r) is the condensate wave function and the second term is the field operator of 
uncondensed particles. The latter is a Bose operator, with the standard commutation relations 

Vi(r),Vl(r)] =5(r-i>). 

It is important to remember that the Fock space J-'i'ip), generated by the operator if), is orthogonal 
to the Fock space .F(if)i), generated by the operator ip u so that after the Bogolubov shift (48) it 
is necessary to work in the space ^(ipi). Mathematical details can be found in Ref. [65]. 

Similarly to property (22), the condensate wave function is orthogonal to the field operator 
of uncondensed particles: 

77*(r)V>i(r) dr = . (49) 



The quantum- number conservation condition, analogous to Eqs. (24) and (43), takes the 
form 

(^(r)) = 0. (50) 

Then Eq. (48) yields 

$(r))=77(r), (51) 
which shows that the condensate function plays the role of an order parameter. 
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The condensate function is normalized to the number of condensed particles 



N = j \v(r)\ 2 dr. (52) 

The number of uncondensed particles gives another normalization condition 

m = (Ni) , (53) 
where the number operator Ni is as in Eq. (26). The total number operator 

N = J ^(r)^(r) dr = N + A\ (54) 

defines the total number of particles 

N = (N) = N + iVx . (55) 

In the Bogolubov representation of the field operator (48), the condensate function and the 
field operator of uncondensed particles are two independent variables, orthogonal to each other. 

4.3 Grand Hamiltonian 

The general self-consistent theory to be presented in this and in the following sections, is based 
on Refs. [63-71], where all details can be found. 

In order to define a representative ensemble, one has to keep in mind the normalization 
conditions (52) and (53). The quantum-number conservation condition (50) is another restriction 
that is necessary to take into account. The latter equation can be rewritten in the standard form 
of a statistical condition by introducing the operator 

A = J [A(r)V>{» + A*(r)^(r)] dr , (56) 

in which A(r) is a complex function that accomplishes the role of a Lagrange multiplier guarantee- 
ing the validity of the conservation condition (50). For this purpose, it is sufficient [71] to choose 
A(r) such that to kill in the Grand Hamiltonian the terms linear in ipi(r). The conservation 
condition (50) can be represented as 

(A) = . (57) 

Taking into account the given statistical conditions (52), (53), and (57) prescribes the form 
of the grand Hamiltonian 

H = H-ti N -HiN 1 -A, (58) 

in which /i and ji\ are the related Lagrange multipliers and H = H[r), tpi] is the energy operator. 
The multiplier /x has the meaning of the condensate chemical potential and /ii can be called the 
chemical potential of uncondensed particles. 

The Hamiltonian average can be represented as 

(H) = (H)-fiN, (59) 

with n being the system chemical potential. Then, from Eq. (58), it follows that the chemical 
potential is 

li = fi n + ii x n x , (60) 
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where the fractions of condensed and uncondensed particles, 

n ° = W ' ni = IT ' (61) 

are introduced. 

It is necessary to stress that the number of Lagrange multipliers in the grand Hamiltonian has 
to be equal to the number of imposed statistical conditions. Only then the statistical ensemble 
will be representative. In the other case, the system would not be uniquely defined. Here, there 
are three conditions, the normalization conditions (52) and (53) and the conservation condition 
(57). 

It is easy to show that the multipliers /io and fi\ do not need to coincide. To this end, let us 
consider the thermodynamic stability condition requiring the extremization of the system free 
energy F = F(T, V, N , Nj, that is, 5F = 0. This gives 

SF = m;, SN ° + w; mi = - (62) 

Substituting here 

"° = aiv ^ = m (63) 

transforms Eq. (62) to the equation 

fi 5N + mSN! = . (64) 

The total number of particles N = N + Ni is assumed to be fixed, so that 5N = and 
5N = -5iVi. Then Eq. (64) reduces to the relation 

(/i - a*i)<Wi = . (65) 

If Ni were arbitrary, then one would have the equivalence of the multipliers /io and fii. However, 
the number of uncondensed particles Ni is fixed for each fixed T, V, and N. That is, SNi = 
and Eq. (65) is satisfied for any multipliers. Hence the multipliers /i and jii do not have to be 
equal. 

It would be possible to say that A^ is fixed for each given T, V, N. But, clearly, this is the 
same as to say that Ni is fixed. In any case, there always exist two normalization conditions 
requiring to introduce two related Lagrange multipliers. 

The Hamiltonian energy operator is 

H = J ^(r) (y-^ + U^j ^) dr + 

+ I J ${r)${i?)<b(T - r')^(r')^(r) drdr' , (66) 

where $(— r) = $(r) is a pair interaction potential and U = U(r,t) is an external potential that, 
generally, can depend on time t. 

Substituting into the grand Hamiltonian (58) the shifted operator (48) results in the form 

4 

h = H(n) . ( 67 ) 

n=0 
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whose terms are classified according to the order of the products of the field operators The 
zero-order term 

if (°) = J r)*(r) (- ^- + U - ^ r)(r) dr + 

+ i y$(r-r / )|^(r / )| 2 |^(r)| 2 rfrrfr / (68) 
does not contain the operators ipi. The first-order term 

H {1) = (69) 
is zero because of the conservation condition (57). The second-order term is 

fi (2) = J ^l(r) (- ^ + U - (i^j ^i(r) dr + 

y $(r-r') [177(^1^1(0^1(0 +r)*(r)r)(r')4(r')Mr) + 

drdr' . (70) 

z, 

Respectively, one has the third-order term 

if ( 3 ) = y$(r-r') ?7*(r)V' 1 t (r / )V'i(r / )V'i(r)+V'i t (r)V'i(r / )V'i(r / )^( r ) rfrrfr ' (71) 



+ 

+ I r ? *(r)^(r')Vi(r / )V'i(r) + ± ^(rMOV'KOV'Kr) 



and the fourth-order term 

H {4) = \ j ^(r)^ 1 t (r / )$(r - r>i(r>i(r) drdr' . (72) 

4.4 Variational Principle 

In the Heisenberg representation, field operators satisfy the Heisenberg equation involving a 
commutator of the operator with the system Hamiltonian. At the same time, in quantum field 
theory, one usually gets the equations for the field operators by extremizing an action functional 
[72,73], which reduces to the variation of the Hamiltonian. Conditions, when these two methods 
are equivalent, are clarified in the following propositions. 

Theorem 6. Let afield operator ip(r) be either Bose or Fermi operator satisfying, respectively, 
the commutation or anticommutation relations 

^(r),^(0] T = <Hr - , Mr), V>(0W = , (73) 

with the upper sign index being for Bose statistics while the lower, for Fermi statistics. Then for 
the products 

m n 

P mn = P+P n , ^ = fl^t( r .) > P n = Y[rKrQ, (74) 

i=i i=i 
where m and n are real integers, one has the commutators 

SP 

Mr), Pmn] = + [(±l) m+n - 1] P m >(r) . (75) 
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Proof. Using the variational derivative 



<ty>t(r) 



5(r - Ti) , 



it is straightforward to find 
5P+ 



^t( r ) 



o(r-ri) , 



5^t( r ) 



5(r - ri)^ t (r 2 ) ± 5(r - r 2 )^ t (ri) , 



5P+ 



5(r - r 1 )^ t (r 2 )V t (r 3 ) ± 5(r - r 2 )^ t (r 1 )V t (r 3 ) + 8(t - r 3 )^ t (r 1 )^(r 2 ) , 



#t( r ) 

and so on. By induction, it follows that 



5P+ 

m 

#t( r ) 



Using the commutator 

[V(rU f (r')] = *(r - r') + (±1 - l)^( r >(r) , 

we derive 

ty(r), P+] = 5(r - n) + (±1 - l)P^(r) , 
ty(r), P 2 +] = 5(r - rO^ra) ± 5(r - r 2 )^( ri ) , 
[V(r), P+] = 5(r - r 1 )^t(r 2 )^t( r3 ) ± 5(r _ r 2 )^( ri )^(r 3 )+ 
+5(r - r 3 )V t (r 1 )V t (r 2 ) + (±1 - l)P+^(r) , 
and so on. From here, using Eq. (76), by induction, we get 

fiP + 

Also, it is easy to check that 

[^(r),P B ] = [(±l) B -l]P B ^(r). 
Then, taking into account that, for any three operators A, B, C, the equality 



A,BC = A,B 



C + B 



A,C 







is valid, we have 

Mr), P mn ] = ty(r), P+] P n + P+[^(r), P n ] . 

Substituting here Eqs. (77) and (78) gives the required Eq. (75). 

Theorem 7. Let F[P mn ] be a linear functional of the products defined in Eq. (74). And 
the linear combination 

f ^ c mn p[p mra ] (i 
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contain only such junctionals for which, in the case of Bose statistics, m and n are arbitrary 
while, for the case of Fermi statistics, m + n is even. Then 



V>(r),F 



5_F 
<^t( r ) 



(80) 



Proof. The proof is straightforward, following immediately from Eq. (75). 

The latter theorem shows that for a large class of functionals the commutator with the field 
operator is equivalent to the variational derivative. The operators of observable quantities are 
in this class, as well as Hamiltonians. This is because for Fermi systems, the field operators 
enter the observables always in pairs, which is necessary for spin conservation. This is why the 
Heisenberg equations for the field operators can be written in two equivalent ways, in the form of 
a commutator, as in the left-hand side of Eq. (80), or in the form of a variational derivative, as 
in the right-hand side of that equation. Note that the standard form of many phenomenological 
evolution equations also involves variational derivatives [74,75]. 

4.5 Evolution Equations 

With the grand Hamiltonian (58), the evolution equations for the field variables rj and ipi read 
as 

• d , , 5H 

I s" (r ' t) = M^y (81) 

for the condensate function, and as 

. d , , , 5H . . 

'a Wr ' iH #lM' (82) 

for the field operator of uncondensed particles. Recall that, in view of Theorem 7, 

5H 



= Mr,t),H]. 



#IM) 

Invoking expression (67) of the grand Hamiltonian gives the equation 

4.(r,t)=(-^ + f /-^),(r,t) + 

+ J $(r - r') [x„(r, r') + X(r, r')] dr' , (83) 

in which the notations are introduced: 

X o (r,r0^(r>?(r>7(r), 

X(r,r') = 4(r')Mr'Hr) + 4(r , Hr , )Mr) + 

+ V *(r')Mr')Mr) + V'lWiWiW ■ (84) 
In these expressions, for brevity, the explicit dependence on time is not shown. 
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Equation (82) yields the equation for the field operator of uncondensed particles: 

+ J $(r - r') [Xx(r, r') + X(r, r')] dr' , (85) 

where 

Xx(ry) ^^(^(rOViW+^^i^W+^KrXrO^W . (86) 

An equation for the condensate function follows from averaging Eq. (81), with the standard 
notation for a statistical average of an operator A as 

(A(t)) = Trp(0)i(t) , 

where p(0) is the statistical operator at the initial time t — 0. So that the condensate-function 
equation is 

• d . , / 5H \ 

Averaging the right-hand side of Eq. (83), we shall need the notations for the single-particle 
density matrix 

pxCr.i/) = (^V)lMr)) (88) 

and the anomalous density matrix 

^(r.r'jE^^W). (89) 
The density of condensed particles is 

p (r) = |77(r)| 2 , (90) 
while the density of uncondensed particles is 

p 1 (r) = p 1 (r,r) = (^ 1 t (r)^ 1 (r)). (91) 
The diagonal element of the anomalous density matrix, 

(7 1 (r) = < 7 1 (r,r) = ^ 1 (r)^ 1 (r)>, (92) 

defines the density of pair-correlated particles as |<7i(r)|. The total density of particles in the 
system is the sum 

p(r)=p (r)+pi(r). (93) 
Also, we shall need the notation for the anomalous triple correlator 

e(r,r') = (^(r')^(r')^(r)). (94) 

Employing these notations gives 

X (r,r')=Po(r>?(r), 
(X(r,r')> =p 1 (r'Mr) + ft (r 1 r'Mr') + ^i(r,r'^(r')+e(r 1 r') ■ 
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Finally, Eqs. (83) and (87) result in the equation for the condensate function 

i l" (r ' t) = (-S + c/ -"»)" (r ' t) + 

+ J $(r — r') [p(r')^(r) + pi(r, r')r)(r') + o"i(r, r')r]*(r') + £(r, r')] dr' . (95) 

Equations for the densities can be obtained from the above equations, with introducing the 
condensate density of current 

jo(M) = " 2^ t^WV^r) -^(r)V^(r)] (96) 
and the current density of uncondensed particles 

JiM) = - ^ (^(r)V^i(r) - [V^j(r)] ^(r)) . (97) 
And let us also define the source term 

T(r,t) = i I $(r-r')[S*(r,r')-S(r,rO]dr', (98) 

with the anomalous correlation function 

S(r, v') = r]*(v) [^(0^(^0+^,0] . 
Then we get the continuity equations for the condensate, 

^Po(r,t) + V-j (r,t) = T(r,t) , (99) 



and for uncondensed particles, 

d 



;;t p 1 (r,t) + V-j 1 (r,t) = -r(r,t). (100) 
The total density (93) satisfies the continuity equation 

^p(r,f) + V.j(r,t) = 0, (101) 

with the total density of current 

j(r,t)=j (r,t)+j 1 (r,t). (102) 
For the anomalous diagonal average (92), we find the equation 
. d 



l dt 



<n(r, t) = 2K(r, t) + 2(U - /^(r, t) + 2 J $(r - r')S(r, r', t) dr' , (103) 
where the average anomalous kinetic-energy density is defined as 

and where we use the notation 

S(r, r, i) = 77(r)77(r>i(r, r) + r]*(r / )^(r) ( j 1 (r, r') + r ] *(r')r ] (r , )a 1 (r) + 
+r ? (r)e(r,r / )+r ? (r / )(^ 1 t (r>i(r)^ 1 (r)) + 
+ ^(^(^(rO^W^W) + (^(r')^(r')^(r)^(r)) + [^(r) + ^(r)]^ - r') . (105) 
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5 Superfluidity in Quantum Systems 
5.1 Superfluid Fraction 

One of the most important features of Bose-condensed systems is superfluidity. Therefore it 
is necessary to have a general definition for calculating the superfluid fraction. Probably, the 
most general such a definition is by identifying the superfluid fraction as the fraction of particles 
nontrivially responding to a velocity boost. 

The systems Hamiltonian H = H[ip] is a functional of the field operator The operator of 
momentum is 

P = J ft(r)inj>(r) dr , (106) 

where p = — iV . 

Boosting the system with a velocity v leads to the Galilean transformation of the field oper- 
ators in the laboratory frame 

i> v (r, t) = ijj(r — vt) exp < i ( mv • r — t ) > , (107) 



expressed through the field operators ip in the frame accompanying the moving system. Then 
the operator of momentum in the frame at rest, 

P v = J ^(r)p^(r) dr , (108) 

P v = | ^(r)(p + mv)^(r) dr = P + mvN . (109) 



transforms into 



Since 



(p + mv) 2 p 2 ///(- 



.2 ™„,2 



= 7T- + V-P + 



2m 2m 2 ' 

the Hamiltonian H v = H[ip v ] for the moving system becomes 

z mv 2 

1 1 1 



//, // • / #(r) ( v • p + — J V(r) dr . (] 10) 

The generalized superfluid fraction is defined through the ratio 

n s (v) - ? { V) " , (HI) 

dv ' ™ v ) 

I V 

in which the statistical averages <•••>„ are determined for the moving system with the Hamil- 
tonian H v , given in Eq. (110). This definition is valid for any system, including nonequilibrium 
and nonuniform systems of arbitrary statistics. 

One usually defines the superfluid fraction for a system at rest, which gives 

n s = limn s (v) . (112) 
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For equilibrium systems, the statistical averages are given by the expressions 



for the moving system, and by 



(A) =^tS (113) 
{A)v = Trexp(-/^) ' (U3) 



for the system at rest. 

In the case of equilibrium systems, the derivatives over parameters can be calculated according 
to the formulas of Ref. [76] . Thus, we have 

= (|;.P„) r /3cov(p^) , (115) 

where the covariance of any two operators, A and B, is 

cov (A, B) =\(AB + BA) v - (A} V (B} V . 

lYiom. Eqs. (109) and (110), one has 

— -P = 3mN ^L=P 
ov av 



Consequently, fraction (111) becomes 



A 2 (P.) 



n a (v) = 1 , (116) 



where the notation for an operator dispersion 

A 2 (A v ) = (A 2 v ) v -(A v )l 
is used. Therefore, for fraction (112), Eq. (116) yields 



n s = l - (117) 
3mNT ' v ; 

with the dispersion given as 

A 2 (i) = (A 2 ) - (A) 2 . 

The quantity 

q = f±A£Z (ii 8 ) 

^ 2mN V ; 

describes the heat dissipated in the considered quantum system. While the dissipated heat in 
the classical case reads as 

Qo = \t. (119) 
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Hence, the superfluid fraction (117) can be represented by the expression 



n s = 1 - £ . (120) 



For an immovable system, the average momentum < P > is zero. Then 

A 2 (P) = (P 2 > (<P> = 0) . 
And the dissipated heat reduces to 



(P 2 ) 

Q = (121) 
^ 2mN V ; 



5.2 Moment of Inertia 



Another way of defining the superfluid fraction is through the system response to rotation. The 
latter is connected with the angular momentum operator 



L 



J V> f (r)(r x p)^(r) dr . (122) 

When the system is rotated with the angular velocity u, the related linear velocity is 

v w = u x r . (123) 
Then, in the laboratory frame, the angular momentum operator takes the form 

K = J V> f (r) [r x (p + mv u )] ij>(r) dr . (124) 

This, using the equality 

r x (w x r) = r 2 cu — (Co ■ r)r , 

gives 

L w = L + m J ${r) [r 2 u - (cu ■ r)r] ^(r) dr . (125) 
The energy Hamiltonian of an immovable system can be written as the sum 

H = K + V (126) 

of the kinetic energy operator 

K = J ft(r) ^- i){r) dr (127) 

and the potential energy part V, respectively. 

Under rotation, the potential energy part does not change, but only the kinetic part varies, 
so that the energy Hamiltonian of a rotating system, in the laboratory frame, becomes 

H u = K + V , (128) 
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with the same potential energy operator V. The kinetic energy operator, in the laboratory frame, 
can be represented [77,78] by the formula 



K U = J ^(r) K *^m UJ ^ (r) * • (129) 



In the rotating frame, where the system is at rest, the kinetic energy operator can be obtained 
from Eq. (129) with replacing u by — uj and, respectively, replacing v w by — v w . Using the 
relations 

(uj x r) 2 = co 2 r 2 — (uj ■ r) 2 , (u x r) ■ p = w ■ (r x p) 
allows us to represent the kinetic energy operator (129) as 

K U = K + + J j V>(r) [u 2 r 2 - (a; • r) 2 ] ^(r) dr . (130) 

Thus, the energy Hamiltonian (128), in the laboratory frame, takes the form 

H uj = H + uj-L + - ^(r) [cA 2 - (w • r) 2 ] ^(r) rfr . (131) 

Rotating systems are characterized by the inertia tensor 

** s a£ (132) 

that, in view of Eq. (125), reads as 

Lf3 = m J ft(r) (r 2 5 aP - r a r p ) ij>(r) dr . (133) 

If one chooses the axis z in the direction of the angular velocity, so that 

uj = ue z , (134) 
then the angular momentum (125) is given by the expression 

L' u = L g + uI gg , (135) 

with the inertia tensor 

L = mJ ^(r) (x 2 + y 2 ) ^(r) dr , (136) 

where the relation r 2 — z 2 = x 2 + y 2 is used. The energy Hamiltonian (128), characterizing the 
system energy in the laboratory frame, can be represented as 

H u = H + uL z + ^i zz , (137) 

with H from Eq. (126). 

The generalized superfluid fraction is defined as 

n» = M^y- • (138) 
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For an equilibrium system, we can again employ the formulas of differentiation over parame- 
ters [76], leading to the derivative 

^<i->„=(^ii)^cov(ii,^| . (139) 
Substituting here 

dK_f dH w _ u 

du 22 ' du u ' 

we come to the expression 

A 2 ( f z ) 

n.( W ) = f - -±-f . (140) 
Considering the superfluid fraction in the nonrotating limit, 



n s = limn s (w) , (141) 



and using the notation 

I zz = Vm{izz) u = mj (x 2 + y 2 ) p(r) dr , (142) 
we obtain the superfluid fraction in the form 

n, = 1 - ^> . (143) 

1 *-zz 

The dispersion of L z is calculated with the Hamiltonian for a nonrotating system. 
Introducing the notation 

I eff = /3A 2 (4) (144) 
allows us to represent the superfluid fraction (143) as 

n s = 1 - ^ . (145) 

For a nonrotating system, one has 

A 2 (L Z ) = (L 2 Z ) ((L z ) = 0). 

Hence I eff = (3 < L 2 Z >. 

5.3 Equivalence of Definitions 

The definitions of the superfluid fraction, considered in Sec. 5.1 and in Sec. 5.2, are equivalent 
with each other. To show this, one can take a cylindrical annulus of radius R, width 5, and 
length L, such that 5 <C R. The volume of this annulus is V ~ 2nRL5. Then the classical inertia 
tensor (142) is I zz ~ mNR 2 . The angular momentum (122) can be written as 

^ f (r) (-i ^ V>(r) dr , (146) 
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where ip is the angle of the cylindrical system of coordinates. 

For the annulus of large radius R, making the round along the annulus circumference, one 
has the path element 51 = R5p. Therefore the angular momentum (146) can be represented as 

L z = RP l , (147) 

being proportional to the momentum 

P^J ft(r) (-i |) ^(r) dr . (148) 
Then the superfluid fraction (143) becomes 

mNT v ; 

The same formula follows from the consideration of Sec. 5.1, if one takes the velocity boost along 
the annulus circumference. 



5.4 Local Superfluidity 

In some cases, it is important to know the spatial distribution of the superfluid fraction that 
would be given by the spatial dependence n s (r). This can be necessary, when one considers 
equilibrium nonuniform systems or systems in local equilibrium [79,80]. 
To describe local superfluidity, we can consider the momentum density 

P(r) = ^(r)p^(r) . (150) 

Following Sec. 5.1, we introduce a velocity boost, which leads to the momentum density 

P ( (r) = ^(r)(p + fnv)^(r) (151) 

in the laboratory frame. The local superfluid fraction is defined as 



n s (r) = lim f v V " v . (152) 



Because of form (151), one has 



' (P w (r) ) B 



— • P„(r) = 3171^ (r)^(r) 



Then the local superfluid fraction (152) reduces to 

.TV 

~3mp{r)T 



„, (r) = ! _ ^mh . m 



Owing to the relation 

p 8 (r) = n s (r)p(r) , (154) 

we get the local superfluid density 

P-M = M- SrnT • (155) 
Integrating the above equation over r and considering the average fraction 

would bring us back to formula (117). 
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5.5 Superfluidity and Condensation 

Usually, Bose-Einstein condensation is accompanied by superfluidity However, there is no 
straightforward relation between these phenomena and the related fractions [3,12]. Thus, in 
two-dimensional systems at finite temperature, there is no Bose condensation, but there can 
exist superfluidity. And in spatially random systems, there can happen local Bose condensation 
without superfluidity. 

The relation between Bose condensation and superfluidity depends on the type of the effective 
particle spectrum and system dimensionality. To illustrate this, let us consider a d-dimensional 
Bose gas with an effective particle spectrum 

cu k = Ak n -fi, (156) 

where A and n are positive parameters and k is <i-dimensional momentum. For the rf-dimensional 
case, the superfluid fraction (117) takes the form 

»■ = 1 - mh ■ (157) 

The integration over the <i-dimensional momenta involves the relation 

dk 2k d ~ 1 dk 



(2%) d (47r) d / 2 r(d/2) ' 



in which T(x) is the gamma function. 

For the condensation temperature, we find 



A 



(An) d l 2 T(d/2)np 



n/d 



(158) 



2T(d/n)((d/n) 

where ((x) is the Riemann zeta function. The latter, can be represented in several forms: 

00 i i roc j.x—1 

c(*) = E - = ttt / 1—r dt ' (159) 

f— ' f Y{x) J e { - 1 



3 



when Rex > 1, and 



i r°° u*' 1 

C(x) = 7 : wT7 — \" / du , (160) 

if Rex > 0. 

Taking into account that T(x) > for x > and ((x) < in the interval < x < 1 tells us 
that there is no condensation for d < n. When d = n, then T c = 0. And T c > for d > n. 
For d > n, the condensate fraction below T c is given by the expression 

/ rp \ d/n 

no = 1 - =- (T<T C ), (161) 



T I 

while the superfluid fraction, under /i — 0, is 



1 _ B( ^ + 2 ~ n ^ T (d+2-n)/n ^ 



n s = l-B(( ^ T ( d+ 2-n)/n ^ (162) 
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where 

2(d + 2 -n)T t x 

~ (4vr) d / 2 r (|) A( rf + 2 )/"mpn 2 rf ' 1 J 

If there is no condensate, then \i is defined by the equation 

2 r W ^ /w( ,)T^ 



in which z = is fugacity and 



1 r zun ~ i A 

9n{z) = ^t-t / — du 



r(n) Jo e" — z 

is the Bose function. The superfluid fraction, in the absence of condensate, is 

1 - Bg {d+2 _ n)/n (z)T( d+2 - n » n . (165) 



n. 



Generally speaking, Bose-Einstein condensation is neither necessary nor sufficient for super- 
fluidity. These phenomena are connected with different system features. Bose condensation 
implies the appearance of coherence in the system, while superfluidity is related to the presence 
of sufficiently strong pair correlations. Thus, there can occur four possibilities, depending on the 
values of the condensate and superfluid fractions: 
(i) incoherent normal fluid 

n = , n s = 0; 



(ii) coherent normal fluid 

(iii) incoherent superfluid 

(iv) coherent superfluid 

n > , n s > 0. 

In this classification, we do not take into account that the system can form a solid [12]. 



no > , n s = 
uq = , n s > 



6 Equilibrium Uniform Systems 
6.1 Information Functional 

The definition of statistical averages involves the use of a statistical operator. The form of 
the latter, in the case of an equilibrium system, can be found from the principle of minimal 
information. This principle requires that, composing an information functional, one has to take 
into account all conditions and constraints that uniquely define the considered system [81]. Only 
then the corresponding statistical ensemble will be representative and will correctly describe 
the system. In the other case, if not all necessary constraints have been taken into account, 
so that the system is not uniquely described, the ensemble is not representative and cannot 
correctly characterize the system. In such a case, one confronts different problems, for instance, 
the occurrence of thermodynamic instability or nonequivalence of ensembles. However, all those 
problems are caused by the use of nonrepresentative ensembles and have nothing to do with 
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physics. A detailed discussion of these problems can be found in Ref. [63]. The construction of 
representative ensembles for Bose-condensed systems is given in Refs. [64,71]. 

A statistical operator p of an equilibrium system should be the minimizer of the Shannon 
information p In p, under given statistical conditions. The first evident condition is the normal- 
ization 

(i) = Trp = 1 , (166) 
with 1 being the unity operator. Then, one defines the internal energy E through the average 

(H) ee Tr P H = E . (167) 

The normalization condition (52) for the condensate function can also be presented in the stan- 
dard form of a statistical condition as 

(N Q ) ee Tr/3A>o = N , (168) 

where N = N 1. Normalization (53), for the number of uncondensed particles, can be written 

as 

(Nx) ee TrpAq = JVi . (169) 
Finally, the conservation condition (57) reads as 

(A) ee TrpA = . (170) 

Note that, in general, the conditional operators do not need to be necessarily commutative 
with the energy operator [80]. For instance, here the operator N commutes with H, but A does 
not have to commute with the latter. 

It is also worth stressing that the average quantities, involved in the statistical conditions, 
do not need to be directly prescribed, but they have to be uniquely defined by fixing other 
thermodynamic parameters. Thus, internal energy is not prescribed directly in either canonical 
or grand canonical ensembles, but it is uniquely defined through the fixed temperature, the 
number of particles in the system, and volume. Similarly, the number of condensed particles 
may be not directly given, but it is uniquely defined, and can be measured, by fixing other 
thermodynamic parameters, temperature, total number of particles, and volume. For confined 
systems, instead of volume, the external potential is given. 

The information functional, under the above conditions, takes the form 

I[p] = Trp hip + A (Trp - 1) + /3(TrpH - E)- 

- /3MTrpA> - N Q ) - Pp^TrpN, - Aq) - /3TrpA , (171) 

in which the corresponding Lagrange multipliers are introduced. Minimizing this functional with 
respect to p yields the statistical operator 

e -m 

?=TZ=SS- (172) 



with the same grand Hamiltonian (58). 
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6.2 Momentum Representation 

For a uniform system, it is convenient to pass to the momentum representation by means of the 
Fourier transformation with plane waves. This is because the plane waves are the natural orbitals 
for a uniform system, which implies that they are the eigenfunctions of the density matrix in the 
sense of eigenproblem (5). 

The field operator of uncondensed particles transforms as 

Mr) = E^^'" > °* = Hff I V>i(r)e- ik ' r dr . (173) 

W k ^ vv J 

We assume that the pair interaction potential is Fourier transformable, 

*(r) = ^ $fce * k ' r > $ * = / H^e'^dr ■ (174) 

k 

The condensate function 77 (r) for a uniform system, is a constant rj, such that 

p (r) = M 2 = po • (175) 

These transformations are substituted into the grand Hamiltonian (67). Then the zero-order 
term (68) becomes 

H<®= Q Po <f - Po) N . (176) 
The first-order term Hi is automatically zero, as in Eq. (69). The second-order term (70) reads 

+ p ($o + $fe) - Pi 
+ \ Yl Po ® k ( a k a -k + a -fc°fc) • (177) 

k^O 

The third-order term (71) is 

^ (3) = VV ^ $p ( a fe° fc +P a -P + a -p a k+P a *) > ( 178 ) 

where in the sum 

k^O, p^O, k + p^0. 
The fourth-order term (72) takes the form 

H(4) = ^Y1 E^M^***-* > ( 179 ) 

q kp 

where 

k^O, p^O, p + q^0, k-q^0. 



as 

r P 

# (2) = E 



k^O 



2m 



a\a k + 
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6.3 Condensate Function 

In the case of an equilibrium system, the condensate function does not depend on time, 

|^,f) = o. 

Therefore, Eq. (95) reduces to the eigenvalue problem 

V 2 



(180) 



2m 



+ U(r) 



rj(r) + 



+ J $(r-0[p(r>(r)+pi(r, r')?7(r') + a^r'Wir') + Z(r,r')]dr' = p rt(r) . (181) 

A uniform system presupposes the absence of a nonuniform external potential. Hence, one 
can set U = 0. The average densities p and p\ are constant. The total particle density is 



p = p( r ) = p + Pl . 



Then Eq. (181) gives 



/i = p$o + / $(r) 



pi(r,0) + cri(r,0) + 



£M) 



182) 



(183) 



The normal density matrix is written as 



1 



Pi( 



( r > r/ ) = T>E n ^ k ' (r ~ r,) 



- /„t. 



0i (r 



ik-(r-r') 



where 

And the anomalous average 
is expressed through 

a k = (a fc a_ fc ) . 

The triple anomalous correlator (94) can be represented as 



n k = (a[a k ) 



e(r,rO = ^E^ k ' (r_r,) 



with 



(184) 

(185) 
(186) 

(187) 
(188) 
(189) 



The diagonal element of Eq. (184) gives the density of uncondensed particles 

Pi = Pi(r,r) = y^rik ■ (190) 
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The diagonal element of the anomalous average (186) is 

a x = <t x (t, r) = ^J2 ak ■ 



(191) 



And the triple correlator (188) leads to 



(192) 



The condensate chemical potential (183) can be rewritten in the form 

Ho = p$o + jj Y] (n k + a k + J • 



(193) 



6.4 Green Functions 



There are several types of Green functions. Here, we shall deal with the causal Green functions 
[81,82] that are called propagators. The set {rj,tj} of the spatial variable rj and time tj will be 
denoted, for brevity, just as j. If there are other internal variables, they can also be included in 
the notation j. 

For a Bose-condensed system, one considers four types of Green functions: 
G n (12) = -i<fy>i(l)*2)> , G 12 (12) = -i(f^(l)Vi(2)) , 

G 21 (12) = -t(f 4(1)4(2)} , G 22 (12) = -t(f4(mm , (194) 
in which T is chronological operator. It is convenient [83] to introduce the retarded interaction 

$(12) = $(n - r 2 )<S(ti - t 2 + 0) . (195) 

Also, one defines the inverse propagators 



G^(12) = 



d V 2 
^ + 2^-^ + ^ 



5(12) - E n (12) , 



^(12) = -E 12 (12) , (7^(12) = -E 21 (12) , 



8(12) - E 22 (12) , 



(196) 



where E aj g(12) is self-energy. Using these, one can write the equations of motion in the matrix 
form 

C 11 1 G'ii + G 12 1 G 2 i = 1 , G 11 1 Gi 2 + G 12 1 G 22 = , 

G 2 \Gu + G 22 G 2 i = , G 2 \G\2 + G22G22 = 1 • (197) 

For a uniform system, when U — 0, one passes to the Fourier transforms of the Green functions 
G a /3(k.,uj), inverse propagators G~p(k, ui), and self-energies E Q/ g(k, w). The inverse propagators 
(196) transform into 

A; 2 

Gn(k,w) = w - h/ii - Sn(k,w) , G7 2 (k,u;) = -Ei 2 (k,u;) , 

2m 
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G , 2 1 1 (k,a;) = -S 21 (k,a;), G 22 l (k, u) = -u - ^ + m - E 22 (k, u) . (198) 

The Green functions enjoy the properties 

G aj8 (-k, u) = G a p(k, u) , Gn(k, -u) = G 22 {k, u) , 

G 12 (k, -w) = G 21 (k,u;) = G 12 (k,a;) . (199) 
And the self-energies also share the same properties 

E a/3 (-k,w) = £ a/3 (k,w) , En(k, -w) = £ 22 (k,w) , 

E 12 (k, -w) = S 21 (k, to) = £ 12 (k, w) . (200) 

Equations (197) yield 

Gll(k,w) - tf(k^) ' Gl2(k ' W) ""^(k^' (201) 

with the denominator 

D(k,u>) = Z 2 12 (k,ou) - G£(k,u)G£(k,u>) . (202) 
6.5 Hugenholtz-Pines Relation 

Hugenholtz and Pines [37], using perturbation theory at zero temperature, found the relation 

//i = En(0,0)-E 12 (0,0) . (203) 

The most general proof of this relation, for any temperature, was given by Bogolubov [16]. 
He proved the theorem, according to which 

|Gn(k,0)|>^, (204) 

where n is the condensate fraction, and 

|G n (k,0)-G 12 (k,0)|>^. (205) 



/,From inequality (204), one has 



lim ]hnD(k,u) = . (206) 



And from inequality (205), it follows that 



2 



^- - /x 1 + E 11 (k,0)-E 12 (k,0) 
2m 



2 



< — ■ (207) 
mn 



The latter inequality leads to the Hugenholtz-Pines relation (203). 

It is important to stress that the expression for fj,i, given by Eq. (203), is exact and, generally, 
it is different from the exact value of /xq in Eq. (183). 
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The Hugenholtz-Pines relation is equivalent to the fact that the particle spectrum is gapless, 
which follows from the following. 

The spectrum e k is given by the zeroes of the Green-function denominator: 

D(k,e k ) = 0, (208) 

which gives the equation 

e k = l - [E n (k,£ fe ) - X 22 (k,e fc )] + _ E? 2 (k, e k ) , (209) 

where 

k 2 1 

u k = — + - Pn(k, e k ) + £ 22 (k, e fc )] - //! . (210) 
2m 2 



In view of condition (206), the limit 



lime fe = (211) 

fc-S>0 



is valid, that is, the spectrum is gapless. 

To find the long-wave spectrum behavior, keeping in mind that the spectrum is uniquely 
defined by Eq. (209), we can use the expansion 



Z a p(k,e k ) ~ £ Q/3 (0,0) + Kpk 2 , (212) 

in which k — > and 

Yf „ = 

Then, defining the sound velocity 



9 



c^^E 12 (0,0) (213) 
and the effective mass 

Tfl 

m * = l + m(S' 11 + S 22 -2S' 12 )' (214) 

we get the acoustic spectrum 

e k ~ck (k ->■ 0) . (215) 

Equation (213), characterizing the general feature of the long-wave spectrum, has been ob- 
tained without approximations, assuming only the validity of expansion (212). Therefore, in a 
Bose-condensed system, the anomalous self-energy S 12 (0, 0) must be nonzero in order to define a 
meaningful nonzero sound velocity. The zero sound velocity would mean the system instability. 
Since expression (213) involves no perturbation theory and no approximations, the condition 

Si 2 (0,0)^0 

is general, as soon as expansion (212) is valid. 
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7 Hartree-Fock-Bogolubov Approximation 
7.1 Nonuniform Matter 

To realize practical calculations, it is necessary to resort to some approximation. The Bogol- 
ubov approximation [13,14] is valid for low temperatures and asymptotically weak interactions. 
The more general approximation, that would be valid for all temperatures and any interaction 
strength, is the Hartree-Fock-Bogolubov (HFB) approximation. Early works [35,36], employing 
this approximation, confronted the inconsistency problem discussed in Sec. 1, because of a gap 
in the particle spectrum. This happened as a result of the use of a nonrepresentative ensemble. 
Employing the representative ensemble of Sec. 4 yields no gap and no any other problems. The 
HFB approximation, applied in the frame of the self-consistent theory of Sec. 4, is gapless and 
conserving [63-71]. 

The HFB approximation simplifies the general Hamiltonian (67). For generality, we consider, 
first, the nonuniform case. 

The third-order term (71) in the HFB approximation is zero. And in the fourth-order term 
(72), the HFB approximation gives 



In what follows, it is convenient to use the notation for the total single-particle density matrix 



4(r)4(r')Mr')Mr) = Pi(r)4(r')M^') + Pi(r>iWi(r) + Pi(r, r)^ 1 t (r> 1 (r) + 

+p 1 (r,r>I(r)Vi(r / ) + ^i(r,r>I(r)^ 1 (r / ) + 
+ ^(r',r)^(r')^(r) - Pl (r) Pl (r') - | Pl (r, r')| 2 - k(r,r')| 2 . 



(216) 



P(r, r') 



77(r)77*(r')+pi(r,r') 



(217) 



and for the total anomalous average 




V(r)r](r') + <r 1 (r,r') . 



(218) 



These equations reduce the grand Hamiltonian (67) to the HFB form 




+ /$(r-r') p(r / )^(r)^ 1 (r)+p(r / ,r)^(r / )^i(r) + 



+ l - <T(r,r')Vl(r')V4(r) + \ a*{v,v')^{v')^{v) dvdv' , 



(219) 



in which the nonoperator term is 



#(0) _ * U( T _ r ') [ Pl (r) Pl (r') + Mim-01 2 + Mr,r')| 2 ] dvdv 1 . (220) 



The condensate-function equation (95) becomes 
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+ I $(r-r')[p(i>(r) + Pl (r, r')^') + ^(r, 0^(0] dr' . (221) 
And the equation of motion (85) for the operator of uncondensed particles now reads as 



■ 9 , . s 



2m 



+ U-/JH) ^i(r) + 



+ J $(r-r') [p(r')^i(r)+p(r, r')^') + a{r,r')4(r')j dr' . 
In the case of an equilibrium system, Eq. (221) reduces to the eigenproblem 

V 2 



(222) 



2m 



+ U r](r) + 



+ J $(r - r') [p(r>(r) + pi(r, r')j)(r) + <7i(r, r')?7*(r')] dr' = /^(r) 
defining the condensate function and the the condensate chemical potential 

1 



(223) 



V 2 

2^ + C/ ' r » 



77 (r) dr + 



+ $(r-r')[po(r)p(r')+pi(r, 077*^)77(0 + ( T 1 (r,r / )r7*(r)7 7 *(r / )]drdr / . (224) 

7.2 Bogolubov Transformations 

The HFB Hamiltonian (219) is a quadratic form with respect to the operators ipi. As any 
quadratic form, it can be diagonalized by means of the Bogolubov canonical transformations, 
whose general properties are described in detail in the book [84] . In the present case, the Bogol- 
ubov transformations read as 

Vi(r) = J2 + t£(r)&t] , h = f [<(r)^(r) - ^(r)^(r)] dv . (225) 

Since ipi is a Bose operator, it should be: 

^K(rK(r')-,*(rK(rl]^(r-r') ) 



K(rK(0 - vl{v)u k {v')] = . 



And, the condition that bk is also a Bose operator leads to the relations 



K( r K( r ) - v t( r ) v p( r )] dr = 6 k P , / [ufc(r)up(r) - u fc (r)iip(r)] rfr = 



(226) 



(227) 



The coefficient functions Uk and are to be defined by the requirement of the Hamiltonian 
diagonalization, under conditions (226) and (227). 
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Let us introduce the notations 
W (r, r') ee 



2m 



+ U(r) -f Jll + J $(r - r>(r') dr' 
+ $(r — r')p(r, r') 



5(r - r') + 



and 



(228) 
(229) 



A(r, r') = $(r - r')<j(r, r') . 
Then the Hamiltonian diagonalization leads to the Bogolubov equations 

J [u(r, r')u k (r') + A(r, v')v k (r')]dr' = e k u k (r) , 

J [u*(r : r')v k (r') + A*(r : r')u k (r')}dr' = -e k v k {v) . 

This is the eigenproblem for the Bogolubov functions u k and v k and the Bogolubov spectrum e k . 
The resulting diagonal Hamiltonian is 



(230) 



with the nonoperator term 



H B = E B + J2 e kb{b k , 



E B = E HFB - £k I \v k (r)\ 2 dr 



(231) 



(232) 



The quasiparticles, described by the operators b k , are called bogolons. Their quantum- number 
distribution is easily calculated, giving 



n k = (blh) = (e fe - 1) 1 , 

which can also be represented as 

n k = - coth ( ^) — 1 . 
2 L \2T ) 

The normal density matrix (88) takes the form 

Pl (r,r') = J2[^u k (r)ul(r') + (1 + jr k )v* k (r)v k (r')} , 

k 

while the anomalous average (89) becomes 

<ri(r, r') = J2^kU k (r)v* k (r') + (1 + n k )v* k (r)u k (r f )] . 

k 

The density of uncondensed particles (91) is 

pi(r) = [*k\u k (r)\ 2 + (1 + n k )\v k (r)\ 2 ] 



(233) 



(234) 



(235) 



(236) 



(237) 
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and the diagonal anomalous average (92) is 

<7i(r) = ^(1 + 2vr fc K(rK(r) . (238) 

k 

The grand thermodynamic potential 

Q = —T In Tre _/3if , (239) 

under Hamiltonian (231), reads as 

fi = E B + T^ln(l-e^ efc ) , (240) 

where the first term, defined in Eq. (232), gives 

Eb = -\J *(r-rO [po(r)p (r')+ 

+ 2p (r)p 1 (r / ) +2? 7 *(r)77(r / )pi(r,r / ) + 2rj* (r)r ) *(r')a 1 (r,r') + 
+ Pl (r)p 1 (r') + |p 1 (r,rO| 2 + | ( 7 1 (r,rO| 2 ]rfrrfr' - J \v k (r)\ 2 dr . (241) 

The above equations are valid for any nonuniform matter, with an arbitrary external potential 
U(v). 

7.3 Uniform Matter 

The previous equations simplify for a uniform case, when there is no external potential. Setting 
U — 0, we can use the Fourier transformation (173) and follow the way of Sec. 6. 
Instead of expressions (228) and (229), we now have 



k 2 



+ p% + p $ fc + 1 n P^k+ P - Hi (242) 



2m V 

p^O 



and 

A fc = po^fc + ^7 Yl ^k+p • (243) 
The HFB Hamiltonian (219) reduces to 

H HFB = E HFB + Y u ka{ak + ^ ^ Ak ( a fc°-fc + a -fe a fc) > ( 244 ) 

with the nonoperator term 

E HFB = H (0) - l - pl%V - 7^^^k+ P (n k n p + a k a p ) , (245) 



2V- 

kp 



in which k^0,p^0. 
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Instead of the Bogolubov canonical transformations (225), one has 

a k = Ukh + v*_ k b ] _ k , h = u* k a k - v* h <J_ k . (246) 
And the Bogolubov equations (230) become 

(u k - e k )u k + A k v k = , A k u k + (u k + e k )v k = . (247) 
The Bogolubov Hamiltonian (231) has the same form, but with 

E B = E HFB + X - J2(e k - u k ) , (248) 

instead of Eq. (232). 

The coefficient functions u k and v k are defined by the Bogolubov equations (247), under 
conditions 

K| 2 - \v- k \ 2 = 1 , u k v* k - v*_ k u_ k = , (249) 

replacing conditions (226) and (227). These functions, due to the system uniformity and isotropy, 
are real and symmetric with respect to the momentum inversion k — > — k. As a result, one has 

2 2 1 2 , 2 A k 

u k -v k = l, u k + v k = — , u k v k = - — , 

2 w k + e k 2 — e k / ot;n x 

% = — • v - = — • (250) 

The Bogolubov spectrum becomes 

e fe = ^ - Ag . (251) 
As is known from Sec. 6, the spectrum has to be gapless, which gives 

= p% + ]- ^2(n k - a k )<5> k . (252) 

k^O 

This is different from the condensate chemical potential (224) that is 

fio = p% + ^ J^(n fc + a k )$ k . (253) 

k^O 

With pi from Eq. (252), expression (242) is 

k 2 1 \ - 

m p^o 

In the long- wave limit, the Bogolubov spectrum (251) is of acoustic form (215), with the sound 
velocity 

c =\M> ( 25s ) 

V m* 
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in which 

A = lim A k = p $ + - V a p % (256) 

p^O 

and the effective mass is 

m * = rr^v 1 ff ^ ' (257) 

where 

$' = — 

p " dp 2 p • 

^From Eqs. (255) and (256), we have 

A = m*c 2 = p $o + ^ Yl ■ ( 258 ) 

Hence, expression (254) can be written as 

k 2 1 v-^ 

Uk = m* c 2 + — + Po($k ~ %) + y 22 n p(®k+ P ~ %) ■ (259) 

Comparing Eqs. (213) and (255) yields 

E 12 (0, 0) = Po $o + 1 ■ ( 26 °) 

And from the Hugenholtz-Pines relation (203), with px from Eq. (252), we get 

£ n (0,0) = (P + Po)$o + ^$^ W (261) 

p^O 

Of course, the same Eqs. (260) and (261) can be derived directly from the Green function 
equations. 

The condensate chemical potential (253) can be written as 

p = S n (0, 0) + E 12 (0, 0) - 2p % ■ (262) 

The difference between Eqs. (252) and (253) takes the form 

p -pi = 2[£ 12 (0,0)-p <l>o] , (263) 

which again tells us that these chemical potentials are different. They coincide only in the 
Bogolubov approximation [13,14], when Ei 2 (0,0) equals po^o- Then p and pi both are also 
equal to po$o and, hence, to each other. 
The momentum distribution (185) is 

(264) 

while the anomalous average (187) reads as 

«* = "t coth (i)' (265) 
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The grand potential (239) enjoys the same form (240), but with 

E B = -^J $(r) [p 2 + 2p oPl (r, 0) + 2p a 1 (r, 0) + 

+ |pi(r,0)| 2 + |<7i(r,0)| 2 ]dr + ^( £fe - ^) , (266) 



k 

which can be transformed to 



E B = ~ y P$o - Po XI ( w p + a p)%- 



7.4 Local-Density Approximation 

When there exists an external potential U(r) and the system is nonuniform, one can use the 
equations from Sec. 7.2. It is also possible to resort to the local-density approximation [1- 
3]. The local-density, or semi-classical, approximation [85,86] is applicable when the external 
potential is sufficiently smooth, such that 



lo dU(r) 



U dr 



< 1 , (268) 



where U and l are the characteristic depth and length of the potential, respectively. 

In this approximation, one looks for the solutions of the Bogolubov equations (230), repre- 
sented as 

u fc (r) = «(k, r) — , v k (r) = v(k, r) — , (269) 

where the functions w(k, r) and i>(k, r) are assumed to be slowly varying as compared to the 
exponentials, so that 

|Vu(k,r)| < fc|u(k,r)| , |Vu(k, r)| < A;|u(k, r)| . (270) 
Then, using the notations 

u(k, r) = |1 + U(r) + 2$ p(r) - ^(r) (271) 
2m 

and 

A^^poW+^rpo, (272) 
one reduces the Bogolubov equations (230) to the form 

[w(k, r) - e(k, r)]u(k, r) + A(r)u(k, r) = , 

A* (r)u(k, r) + [u* (k, r) + e(k, r)]u(k, r) = , (273) 



in which 

<l», 



o = J $(r) dr . (274) 
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The following procedure is analogous to the uniform case. For the coefficient functions, we 
have 



u 2 (k, r) - v 2 (k, r) = 1 , u 2 (k, r) + v 2 (k, r) = 



e(k,r) ' 



u(k, r)v (k, r) = — 



A(r) 
2e(k,r) ' 



u 2 (k,r) 



u(k, r) + e(k, r) 



2e(k,r) 



v 2 (k, r) 



u(k,r) -e(k,r) 
2e(k,r) 



The local Bogolubov spectrum is 



e(k,r) = Va; 2 (k,r)-A 2 (r). 
From the requirement that the spectrum be gapless, 

lime(k, r) = , 

we find 



fc-5>0 



Denoting 

from Eq. (272), we get 
Then Eq. (271) becomes 



A*i(r) = U(r) + [p (r) + 2 Pl (r) - tri(r)]$ 
A(r) = mc 2 (r) , 

mc 2 (r) = [po(r) + <7i(r)]$ • 

A; 2 



w(k, r) = mc (r) 



2m 



The local Bogolubov spectrum (276) takes the form 



e(k,r) = t /c 2 (r)P+^ . 

This shows that c(r) is the local sound velocity. 

With spectrum (282), the bogolon momentum distribution (234) reads as 



?r(k, r) 



coth^^^-1 



In view of the system isotropy, the symmetry properties 

e(-k,r) = e(k,r) , vr(-k, r) = ?r(k, r) 

are valid. 

The single-particle density matrix (235) now transforms into 



ft (r 1 = i^n(k,r)e fc M ) 



(275) 
(276) 

(277) 

(278) 
(279) 
(280) 

(281) 



(282) 



(283) 



(284) 



(285) 
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while the anomalous average (236) becomes 



Here the particle local momentum distribution, replacing Eq. (264), is 



. w(k, r) 
n k, r = — — — r coth 
2e(k,r) 

and, instead of the anomalous average (265), one has 



g(k,r) 
2T 



1 

2 



a(k,r) = - 



mc 2 (r) 



coth 



e(k,r) 



2e(k, r) |_ 2T 
The density of uncondensed particles (237) gives 



and the anomalous average (238) is 



The grand potential (239) reads as 

n = E B + T J In [1 - exp {-/3e(k, r)}] 



rfk 



Here the first term, after the dimensional regularization of the expression 

dk 16m 4 



/ 



e(k, r) - w(k,r)] 



15tt 2 



c 5 (r) , 



takes the form 

Eb = 



(286) 



(287) 



(288) 



(289) 



(290) 



(291) 



(292) 



2 



[ P \y) + 2p (r)p 1 (r) + 2p (r)a 1 (r) + p\{v) + ^(r)] dr + 
8m 4 



+ 



15tt 2 



y c 5 (r)dr 



(293) 



When the system is constrained inside a fixed volume V, then the grand potential f2 = — 
defines the system pressure P = —Q/V, irrespectively of whether the system is uniform or not. 
But, when a nonuniform system is confined inside a trapping potential that does not have rigid 
boundaries constraining the system inside a given volume, then the system pressure cannot be 
defined as —Vt/V. It is possible, being based on the generalized definition of thermodynamic 
limit (13), to introduce an effective volume and effective pressure. However, these quantities are 
different for different potentials and, moreover, they are not uniquely defined even for a given 
potential, hence, they would have no physical meaning. 
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What is well defined for any nonuniform system is the local pressure p(r) that enters the 
grand potential through the equality 

n = - J p(r) dr . (294) 
For the grand potential (291), the local pressure is 

/dk 
In [l-exp{-/Mk, r)}] — + 

+ ^ [p 2 (r) + 2p (r)p 1 (r) + 2p (r)<7 1 (r) + p 2 (r) + ^(r)] - ^ c 5 (r) . (295) 
Equation (295) can be represented as the sum 

p(r) = p (r) +p T (r) , 

in which 

Mr) = [P 2 (r) - rf(r)] i + ^1 _ J> {t) 

and 



Pt(v) = -tJ ln[l - exp{-/3 £ (k,r)}] 



dk 



(2tt) 3 ' 

The latter term, when temperature decreases, tends to zero as 

For asymptotically weak interactions, when $ — >■ 0, Eq. (280), defining the local sound 
velocity, reduces to 

mc 2 (r) ~ p (r)$o • 
In that case, the local pressure (295) simplifies to 

p(r) = i p 2 (r)$ -T J ln[l - exp{-/3 £ (k, r)}] dk 

with the local Bogolubov spectrum 



(2tt) 3 ' 



£(k , r) = ,/ Po ( r) * _ + (_j . 

Such local thermodynamic quantities are common for nonuniform systems, both equilibrium 
[87] and quasiequilibrium [88,89]. 
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7.5 Particle Densities 

In the local-density approximation, it is straightforward to find the densities of particles. Thus, 
the condensate density is 

Po(r) = |^(r)| 2 . (296) 

For an equilibrium system, the condensate function is real. Equation (223) for the condensate 
function, in the local-density approximation, becomes 



V 2 



?7(r) + $ [po(r) + 2 Pl (r) + ^(r)]^!-) = ^(r) . 



(297) 



The simplest way of solving this equation is by means of the Thomas-Fermi approximation, 
when one neglects the spatial derivative, which yields 



Ho-U(r) 

Ptf(i-) = 2pi(r) -<7i(r) . 



(298) 



In the case of cylindrical symmetry, one can introduce the Thomas-Fermi volume Vtf — ttR 2 L, 
with the Thomas-Fermi radius R and the Thomas-Fermi length L defined by the equations 



Ho = U(R, 0) + %[2 Pl (R, 0) + ai{R, 0)] , 



L 



Mo = U [ 0, - j + $ 



2 Pl 0, - + aA 0, 



(299) 



In the Thomas-Fermi approximation, the condensate density is nonzero only inside the Thomas- 
Fermi volume, where 

Po(r)=PTF(r)Q(R-r)Q fL 



— \z 



(300) 



with G(-) being the unit step function. Of course, more correctly, the condensate function should 
be calculated by directly solving Eq. (297). 

The density of uncondensed particles (289) can be written as 



^(k,r) 
e(k,r) 



- 1 



dk 



i 

+ 2 



^(k,r) 
e(k,r) 



coth 



And the anomalous average (290) is 



mc 2 (r) 



coth 



e(k,r) 
2T 



e(k,r) 
2T 



dk 

(2^ 



- 1 



dk 



(301) 



(302) 



At zero temperature, the anomalous average becomes 

If mc 2 (r) dk 
°>V = -2j T^r) 



(303) 



This integral diverges. It can be regularized invoking the dimensional regularization that is well 
defined for asymptotically weak interactions [4]. Employing the dimensional regularization for 
finite interactions requires that the limiting condition 



(7 (r) ^0 (p -)• 0) 



(304) 
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be satisfied [12,66,67,69-71]. This condition takes into account that the anomalous averages and 
Bose condensate always exist together, both being due to the common reason of gauge symmetry 
breaking. As soon as the condensate density is nonzero, the anomalous average is also nonzero. 
And, conversely, when the condensate density becomes zero, the anomalous averages have also 
to disappear. 

Another limiting condition is 



cxo(r) ^ ($ -> 0) . 



(305) 



This condition takes into account that the anomalous average nullifies for the ideal Bose gas 
[12,66,67,69-71]. 

Under conditions (304) and (305), the dimensional regularization gives 



dk 



e(k,r) (2tt) 



^V m$ oPo( r ) . 



7T 



Then Eq. (303) reduces to 



m 2 c 2 (r) 



7T 



V / m$ po(r) • 



(306) 



Thus, at temperatures outside the critical region, the anomalous average (302) can be repre- 
sented in the form 



<7i(r) = (7 - ^ J 



mc 2 (r) 



e(k,r) 



coth 



g(k,r) 
2T 



- 1 



dk 



(307) 



This form can also be used even in the critical region, provided that interactions are weak. 
Strictly speaking, form (307) is valid when one of the following conditions holds true: 



T 

Tr. 



«1, 



«1, 



(308) 



where T c is the critical temperature. 

In the vicinity of the transition point T c , where c(r) — > 0, the anomalous average (302) behaves 

as 

m 2 T 

<7!(r) ~ - — c(r) (T^T C ). (309) 

This behavior guarantees that the Bose condensation transition is of second order for any inter- 
action strength [12,69-71]. 

For the convenience of calculations, the density (301) of uncondensed particles can be trans- 
formed into 



m c (r) 
3tt 2 



2V2 



(Vl + x 2 - l) V2 ^coth 



mc 2 (r) 
2T 



x 



1 dx 



(310) 



and the anomalous average (307), into 



0-1 (r) = oo(r) 



m 3 c 3 (r) 



2V2tt 2 



1) 



1/2 



Vl + x 2 



coth 



mc 2 (r) 
2T 



x 



Ibis. 



(311) 



The sound velocity here is defined by Eq. (280). 
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The local superfluid density has been introduced in Eq. (155). In the local-density approxi- 
mation, for an equilibrium system, we have 

Ps(r) = P (r) - ^ , (312) 

with the local dissipated heat 

— [n(k, r) + n 2 (k, r) - a\k, r)] _ . (313) 

In view of Eqs. (287) and (288), this yields 

1 r k 4 dk 



(4n) 2 mJ sinh 2 [5(k,r)/2T] ' 
which can be transformed into 

= mV(r) (v^+^-l) 3/2 ^ 
1 J (2ir) 2 V2j v / rT^sinh 2 [mc 2 (r)x/2T] ' 1 ' 

It is necessary to stress the importance of taking account of the anomalous average. If in Eq. 
(313), one would omit this anomalous average, then the dissipated heat would be infinite, hence 
the superfluid density would not exist at all. But, taking the anomalous average into account 
renders the dissipated heat (314) a well defined finite quantity. The fact that the anomalous 
average is crucially important for describing superfluidity should be apparent remembering that 
|<Ji(r)| is the density of pair-correlated particles. These pair correlations are, actually, responsible 
for the existence of superfluidity as such. Therefore, when there are no pair correlations, there 
is no supefluidity. 

Having all particle densities defined makes it possible to study their spatial distributions and 
to calculate the average condensate, n , and superfluid, n s , fractions, as well as the fraction n\ 
of uncondensed particles, given by the equations 



Ul N 



^ J pi(r) dr , n + n x = 1 . (315) 



8 Local Interaction Potential 
8.1 Grand Hamiltonian 

Till now, the consideration, for generality, has been accomplished for any type of the symmetric 
interaction potential $(— r) = $(r), with the sole restriction that this potential be integrable, 
such that integral (274), defining $o, be finite. 

When particles interact with each other through a potential, whose effective interaction radius 
r is much shorter than the mean interparticle distance a, then this potential can be represented 
in the local form 

$(r) = $ 5(r) , $ = 4tt ^ , (316) 
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in which the interaction strength <J> is expressed through s-wave scattering length a s and mass 
m. For uniform systems, the potential is called stable [90] when $ is positive. For trapped 
atoms, a finite system can be stable also for negative interactions [1-3,91]. 

The grand Hamiltonian (67), for the local interaction potential (316), contains the following 
terms. The zero-order term (68) reads as 



= yV(r) (- ^ + U - /i ) V(r) dr + yJ |r?(r)|4rfr • 
The first-order term, as always, is zero. The second-order term (70) is 

2|??(r)| 2 v4(rM(r) + ^(r))^)^) + \(v(r)M(r)4(r) 



(317) 



+ $o 

The third-order term (71) becomes 



dr 



H (3) = $ 



?y*(r)^J(r)^i(r)^i(r) + v4( r )^i( r )?M r M r 



dr 



(318) 



(319) 



And the fourth-order term (72) reduces to 

H (4) = y J 1>l(r)4(r)Mr)Mr) dr. 
8.2 Evolution Equations 



(320) 



Evolution equations, derived in Sec. 4.5, simplify for the local potential (316). The same nota- 
tions (90) to (93) can be used. But, instead of (94), we define 



£(r) = (V>I(r)Vi(r)Vi(r)) . 
In addition, we shall employ the notation 

£i(r) ee (^(r)^(r)Vi(r)) . 
Equation (95) for the condensate function yields 



(321) 



(322) 



(323) 



+ $ [po(r)?7(r) + 2p 1 (r)?7(r) + (7 1 (r)?7*(r)+^(r)] . 
The continuity equations (99) to (101) have the same form, but with the source term 

r(r,t)=^ [S*(r)-S(r)] , (324) 

with the anomalous correlation function 

S(r)=^(r)[^(rV 1 (r)+e(r)]. 
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Equation (85) for the operator of uncondensed particles changes to 



■ 9 , . , 



2m 



+ U - m) V>i(r, t) + $o Xi(r, r) + X(r, r) 



(325) 



Equation (103) for the anomalous average becomes 

. d 



i — ^(r, t) = 2K(r) + 2(U - ^(r) 



+ 2$ 



r, 2 {v) Pl {r) + 2p (r)(7 1 (r) + 2??(r)£(r) + ?f (r)a(r) + (^{(r)^i(r)^i(r)^i(r)> 

+ 2[77 2 (r) + <xi(r)]$(0). 



+ 



(326) 



The quantity $(0), under the local potential (316), is not defined and requires to be specified by 
additional constraints. 

A straightforward formal way of giving some meaning to this quantity would be by remem- 
bering that the delta potential (316) is the limiting form of a potential with a finite interaction 
range r , such that r a. For instance, potential (316) could be treated as the limiting form 
of the potential 

/ 3 r 2^ 
$(r) = Aexp - 



where r — > 0, so that the integral 



2rl 



$(r) dr 



(327) 



is fixed as in Eq. (316). The interaction radius is defined as 

rl = — l r 2 $(r) dr . 
$o J 

These requirements give 



(328) 



2tt 



3/2 



Then, for potential (327), the quantity $(0) should be defined as 

$(0) = A 



ir mr, 



3 • 



(329) 



However, we have to always remember that the local interaction potential (316) is an effective 
potential modeling particle interactions for the processes occurring at the interparticle distance 
much larger than the interaction radius. In order to characterize the processes at short distance, 
one has to use a different effective potential that takes into account particle correlations [5,81,83]. 
The latter, in particular, show that two particles cannot exist at the same spatial point. This is 
equivalent to saying that $(0) must be set to zero. 

The four-operator term can be simplified as 



(V4(r)^i(r)^i(r)^i(r)) = 3pi(r)<7i(r) , 
while the three-operator terms are left untouched. 



(330) 
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Then the evolution equation (103) for the anomalous average leads to 



% |- <7!(r) = 2K{v) + 2{U - /^(r) + 
at 

+ 2$ ^ 2 (r) Pl (r) + 2 Po (r) ( 7 1 (r) + Zp^a^r) + 2??(r)£(r) + ^(r)^ (r)] , (331) 

with the anomalous kinetic-energy density K(r) given by Eq. (104) and with $(0) set to zero. 

The derived evolution equations can be used for studying the initiation of Bose-Einstein 
condensation and also the decoherence processes in systems with spontaneous symmetry breaking. 
For finite systems, with N degrees of freedom, coherence can persist [92-95] during the time not 
longer than that of order N/T. 

8.3 Equilibrium Systems 

Considering equilibrium systems, we follow Sees. 6 and 7, substituting in the corresponding 
equations the local potential (316). The grand Hamiltonian (219), in the HFB approximation, 
reads as 

Hhfb = E HFB + j ^(r) (- ^- + U - p^j ^i(r) dr + 



+ $0 

where the nonoperator term is 



2p(r)Vl(r)Vi(r) + \ a(r)4(r)4(r) + l - <7*(r)V>i(r)V>i(r) 



dr , (332) 



Ehfb = # {0) " y | [2p?(r) + a 2 (r)] dr , (333) 

and the notation 

p(r) =p (r)+pi(r) , <r(r) = ^ 2 (r) + ^(r) (334) 
is used. The condensate-function equation (223) becomes 

V 2 \ 

— + UJ ri(r) + $ (Po(r)^(r) + 2p 1 (r)? 7 (r) + <7i(r)V(r) + £(r)] = p r](r) . (335) 
For the condensate chemical potential (224), we have 

+ Wj [P ° W + 2 <" Mpl(r) + (I'M? M*)] dr . (336) 
Employing the Bogolubov transformations (225) yields the Bogolubov equations 

w(r)u fc (r) + A(r)u fc (r) = e k u k (r) , u(r)v k (r) + A*(r)w fc (r) = -e k v k {v) (337) 
replacing Eqs. (230), with the operator 

w(r) = - ^ + C/(r) - ^ + 2$ p(r) , (338) 
2m 
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instead of Eq. (228), and with 

A(r) = $ [»7 2 (r)+<7i(r)] , (339) 

instead of Eq. (229). The chemical potential pi is defined by the requirement that the spectrum 
Ek be gapless. 

The grand potential has the same form (240), with 

E b = ~yJ [Po(r)+4po(r)p 1 (r) + 2p?(r) + 2(^(r))V 1 (r) + ( T 1 2 (r)] dr - 

~ I £ kW{v)\ 2 dr , (340) 
k J 

in agreement with Eq. (241). 
8.4 Uniform Systems 

Resorting to the Fourier transformation (173) gives the following terms of the grand Hamiltonian 
(67). The zero-order term (317) is 

H (0) = (^_ p 2 % _^ p ^y (341) 

The first-order term is, as always, zero. The second-order term (318) reads as 

H {2) = + 2 A>$o - Pi) aW + \ Po$o { a l a -k + a -^ a k) • (342) 

The third-order term (319) yields 

H {3) = M^oY^ ( a W +P a- P + al p a\ +p a k ) , (343) 

' hp 

where 

k^O, k + p^0, p^O. 
And the fourth-order term (320) becomes 

q kp 

in which 

MO, p^O, p + q^0, k-q^0. 
Equation (183) defines the condensate chemical potential 



p = P$o + 77 ( n k + & k + -7= J • 



(345) 
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The grand Hamiltonian (332) in the HFB approximation transforms to 

A 

fc^O k^O 

where the nonoperator term (333) is given by the expression 



Hhfb = Ehfb + ^2 w fc4 a fc + y { a l a -k + > ( 34 6) 



E}JFB $0 / 2 o 2 , 2\ / 0/ ,rA 

— — = y ~ 2/?1 + a i) ~ • ( 34T ) 

Also, here the notation 

k 2 

w fc = A*i + 2p$ (348) 

2m 

is used, and, instead of Eq. (339), we have 

A = $ (po + ai) ■ (349) 

The diagonalization of Hamiltonian (346) is done by the Bogolubov canonical transformations 
(246), resulting in the Hamiltonian of the Bogolubov form (231). 
The condensate chemical potential (253), or (336), reads as 

Po = $o(po + 2pi + <7i) . (350) 

And the chemical potential (252) becomes 

Pi = $ (po + 2pi - <7i) . (351) 

Using the latter in Eq. (348) gives 

k 2 

co k = — + $o(po + ^i)- (352) 
2m 

The long-wave spectrum is acoustic, 

e k ~ck (k ->■ 0) , (353) 

with the sound velocity 

c=\—, (354) 
V m 

which is defined by the equation 

mc 2 = $ (po + <ri) . (355) 
Combining Eqs. (352) and (355) yields 

k 2 

Uk = mc 2 + - — . (356) 
The solution to the Bogolubov equations (247) results in the Bogolubov spectrum 



/ k 2 x 2 



\ (<*) 2 - J (357) 
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Calculating Eq. (340), we resort to the dimensional regularization giving 

dk 16m 4 c 5 



/ 



(e k — oj k ) 



15tt 2 



Then Eq. (340) reduces to 



m 4 c 5 



The system pressure can be expressed through the grand potential (240), which gives 



E 



V 



V 



-fit* 



dk 

(2^p " 



(358) 



(359) 



The integral in Eq. (359) corresponds to thermal pressure. It can be calculated by transforming 
it to the form 



J K ' (2vr) 3 2v/27r 2 Jo 



In 



mc 

1 - exp ( - — x 



(Vl + x 2 — l) xdx 



At low temperatures, such that T <C mc 2 , one can expand the integral as 



J In (1 - e - pe ") 



dk 



[mc) 



2ix 2 



T \ 3 15 / T x ' 



2 \ mc 2 



Therefore, the zero-temperature pressure is 



P=-^ ( T = °)- 
The internal energy is given by the expression 

E=(H)+/iN , 



(360) 



(361) 



in which the average of the grand Hamiltonian is as in Eq. (59) and the system chemical potential 
is defined in Eq. (60). At zero temperature, the internal energy (361) yields the ground-state 
energy 

E = E B + fiN (T = 0) , (362) 

where we take into account that 

(H) = (H B ) = E B (T = 0) . 
Then the ground-state energy is given by the equation 



4^5 



E 2%a s , 2 2 2 x 8m 4 c 

(P + Pi - 2pi<ri - crj + 



iV mp vr ! ' ' u 15tt 2 p 

For convenience, let us introduce the dimensionless ground-state energy 

2mE 



e 



1 ATp3/2 
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(363) 



(364) 



7 = P 1/3 a s 



(365) 



and the dimensionless gas parameter 
The ground-state energy (364) at weak interactions, when 7 <C 1, allows [69,70] for the expansion 

e ~ 47T7 + ^ 7 5/2 + ^ 7 4 • (366) 
15 9 

The first two terms here reproduce the Lee-Huang- Yang result [96-98] . When particle interactions 
are strong, so that 7>1, then [69] one has 



eo .8. 7+ | M ' /3 -5(3^) i/3 i + ^(3,r 3 ^. 



(367) 



8.5 Atomic Fractions 



For the local interaction potential (316), it is straightforward to calculate all atomic densities 
and fractions. The condensate density 



Po = P ~ Pi 

is expressed through the density of uncondensed particles 



2e k 



coth (i*-) - 1 

\2T ) 2 



dk 



that can be rewritten as 



The anomalous average 



-s 



dk 



dk 



+ 



2e k 



coth ( 
\2T 



dk 



(368) 



(369) 



(370) 



(27T) 



mc 



/ e fe \ dk 
2e k C ° th V27V (2^ 



(371) 



can be treated as in Sec. 7, by separating the term 

mc 2 dk 



0o 



2e fc (2tt)3 ' 



which gives 



o\ = a - 



mc 
2^7 



coth 



dk 



(—)-- 

\2TJ J (2tt) 3 



(372) 



(373) 



Expression (372) diverges, but can be regularized invoking the dimensional regularization, as 
in Sec. 7, resulting in 

1 dk 2m I — 

— 7T^3 = ~ V m Po*o • 



7T Z 



Then Eq. (372) becomes 



(374) 
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The anomalous average (373), with the separated term (374), is valid if at least one of conditions 
(308) is satisfied. 

In the vicinity of the condensation temperature T c , it is necessary to use the anomalous 
average in the form 

u\ — <r-ia, 



2tt 

which follows from Eq. (371) and guarantees the second-order phase transition. 
The density of uncondensed particles (370) can be represented as 



Pi = 



m 3 c 3 
3^ 



1 + 57if (^-0 



1/2 



coth 



and the anomalous average (373), as 



<7\ = cr - 



m 3 c 3 r(VTT^-i) 



2V2n 2 Jo 



1/2 



coth 



mc 



2T 



mc 
~2T 



x 



dx 



x - 1 



dx 



The superfluid fraction 



n, = 1 - 



2Q 
3T 



is expressed through the dissipated heat 



Q 



A 2 (P) _ (P 2 ) 
2mN ~ 2mN 



In the HFB approximation, we get 



which reduces to 



Q 



Q 



2^( nfc + % -^(2^' 



1 ( 



k A dk 



sinh 2 (e fe /2T) 



The superfluid fraction (378) leads to the superfluid density 

(Vl + x 2 — l) 3//2 xdx 



[mc) 



Ps = P~ 



6V2~7i 2 mT Jo y/l + x 2 smh 2 (mc 2 x/2T) 
The particle densities are related to particle fractions 



_ Po n _ Pi 

no = — = 1 — ni , rii = — 

P P 
Also, let us define the dimensionless anomalous average 

a = — 



n, = 



_ Ps 



and the dimensionless sound velocity 



mc 



,1/3 • 



(375) 



(376) 



(377) 



(378) 



(379) 



(380) 



(381) 



(382) 
(383) 

(384) 
(385) 
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At zero temperature, we have 



s 3 s 3 



a = ^ V^~o (T = 0) . (386) 

When interactions are weak, such that 7 < 1, then the following expansions are valid for the 
condensate fraction 

8 o /2 64 , 640 0/2 , „v 

^'-Vi^'-S 1 -^^' (387) 

sound velocity 

o r- 1/2 16 2 32 7/2 3904 5 /000 , 

. * 2^ + y T + V? 7 " 27T 1 ' (388) 

and the anomalous average 

- - 4= ^ + ~^-% 7 9/2 • (389) 

As is seen, the anomalous average is three times larger than the normal fraction of uncondensed 
particles: 

— ~ 3 (7 < 1) . (390) 

Hi 

This emphasizes again that the anomalous average in no way can be neglected for a Bose- 
condensed system. 

For strong interactions, when 7 ^> 1, we find the following expansions for the condensate 
fraction 

TT 1 1 fn 5 \ 1/3 1 

n ~ 



64 7 3 512 V 9 J 7 5 ' (391) 
sound velocity 

, 0,1/3 1 /tt 5 \ 1/3 1 1 /tt 7 \ 1/3 1 , x 
and the anomalous average 

(97T) 1 / 3 1 _ ^ J_ _ J_ 1 J_ J_ 

a " 1 7 64 7 3 128 V 3 7 ' .".12 \<) ) 7 5 ' ( ' i9 ' 3) 



Though now the anomalous average is smaller than n 1; but it is much larger than the condensate 
fraction: 

— ~ 15.5167 2 (7 > 1) . (394) 
n 

Thus, the anomalous average is always of crucial importance and can never be neglected. 
At low temperatures, such that 

« 1 , (395) 
we find [69-71] the fraction of uncondensed particles 

(mc) 3 (mc) 3 ( T x> " 



* * t£ + w 1 • < 396 > 
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anomalous average 

(T (mc) 3 / T 
p lip \mc 
condensate fraction 

(mc) 3 (mc) 3 / T x - " 
3tt 2 p Yip \mc 

and the superfTuid fraction 

1 27r 2 (mc) 3 f T s ' 



.. , , (399) 
45p \mc z J 

Notice that the temperature corrections for a are the same as for n\ and Uq. 
Bose-Einstein condensation happens at the temperature 

2/3 

(400) 



nn 27T 

T — 

J- c — 

m 



P 



C(3/2) 



In the critical region, where T — > T c , so that 

< 1 , (401) 



mc 2 



T 

we find [69-71] the expansions 

/T\ 3/2 (mc) 3 m 2 cT /T\ 3/2 (mc) 3 

The superfluid fraction disappears together with the condensate fraction. 

Passing to dimensionless quantities, it is convenient to consider the temperature deviation 

r = l-^r (T<T C ) (403) 

from the dimensionless transition temperature 

mT 27r 



Then we obtain 



3?r 9tt / 2?r ' 
s — — r + i — 1 — — 9 r > ni~ 1— - t + -t 



f c 4t c V 7*c/ 2 

3 3 / 67r \ 2 3 3 2 

°" — ~o r +o 1 + — T^ r ' no ~ - r - - r 



2 8 V ill) 2 8 

3 / 132.411 



c 



-2' 8V 1 + ^ll^» r2 - (405 > 



As is evident, though the anomalous average tends to zero, as T — > T c , but it is of the 
same order as the condensate fraction, hence, again the anomalous average cannot be omitted. 
If it were neglected, the transition would become of first order, which is principally incorrect 
[67,71]. While accurately taking account of the anomalous average renders the Bose-Einstein 
condensation the correct second-order transition, as is obvious from expansions (405). 
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9 Disordered Bose Systems 



9.1 Random Potentials 

The properties of Bose systems can be essentially changed by imposing external spatially random 
potentials. In this section, the theory is presented for the case when such random potentials are 
imposed on a uniform system. There is vast literature studying Bose systems inside randomly 
perturbed periodic lattices (see the review article [12] and recent works [99-101], where further 
references can be found). But this is a different problem that is not touched in the present section. 
In several papers (e.g. [102-106]) the influence of weak disorder on a uniform Bose-condensed 
system has been studied. Strong disorder can be treated by means of numerical Monte Carlo 
simulations [107]. In the present section, the analytical theory is described, which is valid for 
arbitrarily strong disorder. The consideration below is based on Refs. [108-110]. 
The system is described by the grand Hamiltonian 

H — H — fi N - HiNi - A , (406) 

with the energy Hamiltonian 

H = J ft{r) - + U(r) + f(r) ftr) dr + y J ft (r)ft (r)ftr)ftr) dr (407) 

containing a random external potential £(r). Other notations are the same as in the previous 
sections. 

The random potential, without the loss of generality, can be treated as zero-centered, such 
that 

«e(r)» = . (408) 

The double brackets imply the related stochastic averaging [111]. The random-potential correla- 
tions are characterized by the correlation function 

<(£(r)£(r')>) = i?(r-r'). (409) 

One can use the Fourier transformations 

R(r) = ^J2 Rk ^' r > R *=[ R(r)^' rdr ■ ( 41 0) 

k 

Then Eq. (409) yields the correlators 

((C k Q) = S kp R k , = S- kp R k ■ (411) 

The quantum statistical averaging, involving a Hamiltonian H, for an operator A, is denoted 

as 

The total averaging, including both the quantum and stochastic averagings, is denoted as 

(i) = «(<i>*)». (413) 
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The grand thermodynamic potential is given by the expression 

D, = -T((lnTre-^ H )) , (414) 

corresponding to the frozen disorder. 

In addition to the particle densities, considered in the previous sections, for a random system, 
it is necessary to introduce one more density. This is the glassy density [108] 

Pg^\j [((\(Mr)) H \ 2 ))dr. (415) 



With the Fourier transform 



we come to 



V 



^ V k^O 



Pg = ^E«I^! 2 ))' ( 416 ) 



V 

k^O 



where 



a k = (a k ) H ■ (417) 

Because of condition (50), we have 

(R» = (a*) = . (418) 
However, quantity (417) is not zero. The glassy fraction is given by 

nc = P -f = ± J((\{Mr)) H \ 2 ))dv, (419) 

which can be represented as 

"a = ^£«W»- (420) 

To better illustrate the idea of the approach we aim at developing, let us set U = 0. This will 
simplify the consideration. Then the grand Hamiltonian writes as 

4 

H = J^# (n) + Hi. , (421) 

n=0 

where the first sum consists of terms (341) to (344), while the last part 

= Poto^/V + V^Y1 { a ^ k + ^ ak ) + \ Yl a l a p^- P (422) 

k^O fcp (^ ) 

is due to the presence of the random potential. 
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9.2 Stochastic Decoupling 

The sum in Hamiltonian (421) can be treated in the standard way by resorting to the HFB 
approximation, as in the previous sections. But the part (422), characterizing the action on 
particles of the random potential, has to be treated with caution. If one would apply to this part 
the simple HFB-type approximation 

a{ a p£k-p ->■ ( a k a p)£k- P + a \ a p{ik- P ) - (a{ a p) (£k- P ) , 
then the influence of this part, because of Eq. (408), would reduce to the trivial mean-field form 

Yl a k a p£k-p -> Pl£o w , 
fcp(^O) 

containing no nontrivial information on the action of the random potential on particles. 

In order not to loose the information on the influence of the random potential, we employ the 
idea of stochastic decoupling that has been used earlier for taking into account stochastic effects 
in different systems, such as resonant atoms [112-115] and spin assemblies [116-121]. 

In the present case, the idea is that the simplification of the third-order expression in the last 
term of Eq. (422) should include only the quantum statistical averaging, but not the stochastic 
averaging, thus retaining undisturbed stochastic correlations. This idea can be represented in 
several equivalent ways. We can write 

(4 a p£fc- P ) = ({a* k a p £ k _ p }} , (423) 

which is equivalent to 

(4 a p£fc- P )tf = afctt^fc-p • (424) 
In turn, the latter is equivalent to the decoupling 

a\a p = a\a p + a* k a p — a* k a p . (425) 

Then, we introduce [108-110] the nonuniform canonical transformation 

a k = u k b k + v*_ k b ] _ k + w k if k , (426) 

whose coefficient functions are defined so that to diagonalize the grand Hamiltonian (421) in 
terms of the operators b k . The latter are treated as quantum variables with the condition 

(b k ) H = . (427) 

The variables cp k represent stochastic fields. In view of Eqs. (426) and (427), we have 

ttfc = (afc)fr = w k (p k . (428) 

Diagonalizing Hamiltonian (421) results in the relations 

2 _ Wk + gfc 2 _ U k — gfc 

Uk ~ 2e k ' Vk ~ 2e k ' 

771 I 

u k v k = - - — , w k = - ■ , (429) 

2e k 2u k + mc z 
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in which 

k 2 

oo k = h mc 2 . (430) 

2m 

The Bogolubov spectrum 



* = \l (<*)*+ I £i) 2 ( 431 ) 



has the standard form, with the sound velocity defined by the equation 



P 

The stochastic field satisfies the Fredholm equation 



mc 2 = I n + — ) p$ • (432) 



p 



Hamiltonian (421) acquires the diagonal form 

H = E B + J2 zAh + H ran , (434) 
k 

where the last term 

Hran = V?0 \fN~0 (435) 

characterizes the explicit influence of the random potential on the system energy. 
For the particle momentum distribution (185), we get 

Uk = £ C ° th (If) " l + ^ a ^> ( 436 ) 
and for the anomalous average (187), 

= ~ ^ coth (^) + ((K| 2 » . (437) 

Expressions (436) and (437) possess, as compared with Eqs. (264) and (265), additional terms 
caused by the random potential. From Eqs. (428) and (429), it follows that 

«M 2 » = , (<l r' 2> L ■ («8) 
(ou k + mc z y 

The partial chemical potentials (350) and (351) are of the same form 

«o = (p + Pi + , Mi = (p + Pi ~ <7i)$ ■ (439) 

But the quantities 

entering them are now different. The density of uncondensed particles becomes the sum of two 
terms, 

Pi= Pn + Pg ■ (440) 
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The first term is the normal density 

dk 



Pn = - — coth ( ^) — 1 
I N 2 J e k \2TJ 



(2tt) 3 

that, as earlier, is due to finite temperature and interactions, while the second term, 



(441) 



n f ((M 2 )) dk 
J (u k + mc 2 ) 2 (2tt) 6 

is the glassy density produced by the random potential. 
The anomalous average is also the sum of two terms, 

o\ = cr N + p G ■ (443) 

The first term is 

1 f mc 2 , / Su \ dk , . . .s 

°" = -2j — coth iw) ' (444) 

while the second term, caused by the presence of the random potential, coincides with the glassy 
density (442). 

The partial chemical potentials (439), with expressions (440) and (443), become 

A*o = (p + Pn + a N + 2p G )® , fa = (p + p N - <t n )®o , (445) 

essentially differing from each other. 

The superfluid density (378) requires the knowledge of the dissipated heat (379). The latter 
also reduces to the two-term sum 

Q = Qn + Qg ■ (446) 
The first term is analogous to Eq. (380) giving 

1 f k 2 dk 

Qn ~ 8mp2 J sinh 2 (e fc /2T) (2nf ' [ V 

And the second term 

VG 2mp J e k (u k + mc 2 ) \2T ) (2tt) 3 v ; 

is the heat dissipated by the glassy fraction. Thus, the superfluid density (378) takes the form 

9.3 Perturbation-Theory Failure 

Considering the case of weak disorder, it is tempting to resort to perturbation theory with respect 
to disorder strength. In doing this, one has to keep in mind that such a perturbation theory can 
fail. Therefore, the results derived by means of perturbation theory may be not reliable. To 
illustrate this, let us consider the average energy 

Eran = (H ran ) = ( (y9 > ) v^o , (450) 
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related to the random term (435). 

Assuming that disorder is weak, one could think that Eq. (433) could be treated perturba- 
tively, by means of the iteration procedure starting with 

<4 0) = VPo & ■ 

The first iteration gives 

(1) , c jPo \^ 6s-p£p 

V V ui v + mc 2 
V v 

Using this in Eq. (450) yields 

«ig 2 » 



Eil = -Po J2 



V 



u v + mc 2 



In view of correlators (411), one gets 



EV =- [ ^ . (451) 

J cu p + mc 2 (2vr) 3 1 } 

That is, the direct influence of the random potential would lead to the decrease of the system 
energy. It is exactly this expression (451) that has been obtained by several authors (see, e.g., 
[122]) employing perturbation theory. 

However, from Eqs. (418) and (428), involving no perturbation theory, it is seen that 

(K» = o , «^» = o . 

Consequently, the random energy (450) is exactly zero: 

Eran = (H ran ) = . (452) 

Also, using perturbation theory in calculating sound velocity, some authors (e.g. [122]) find 
that the speed of sound would increase due to the random potential. Contrary to this, in our 
theory [108-110], the sound velocity decreases, which looks more natural and is in agreement with 
other [105] calculations. Really, it looks to be clear that the occurrence of an additional random 
potential should lead to additional scattering and, hence, to the decrease of sound velocity. 
Some other contradictions resulting from the use of perturbation theory are illustrated in Refs. 
[108,110]. 

Note that all consideration above is valid for any type of disorder characterized by the corre- 
sponding correlation function (409). 



9.4 Local Correlations 

To proceed further, let us consider local correlations, described by the delta-correlated disorder, 
when the correlation function (409) is 

R(r) = R 5(r) . (453) 

Then Eq. (411) gives 

((e k Q) = h P Ro ■ (454) 
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The solution to the Fredholm equation (433) can be well approximated [108] by 



1 _i_ J_ v ^ 

y/y L^ v uip+mc 2 

It is convenient to introduce [108] the disorder parameter 



= : , ^— • (455) 



hoc p i/6 hoc 

being the ratio of the mean interparticle distance versus the localization length 

l ^^k- (457) 

Employing the self-similar approximation theory [123-132] results in 

«w 2 » = ■ < 458 > 

where s is the dimensionless sound velocity (385). 

The local disorder, with the delta correlation (453), allows for more straightforward calcula- 
tions. At the same time, it gives good understanding of the influence of disorder on the system 
even for the general case of nonlocal disorder. If the random potential £(r) is characterized by a 
finite strength Vr, with the correlation function (409) having a finite correlation length Ir, then, 
to pass to that case, one should make the replacement 

R = V^l\ , (459) 

which follows directly from the definition of correlator (409). As a result, the localization length 
(457) becomes 

Aix 

l ioc = tt-^TTIs ■ (460) 



The particle fractions 



— P° _ Pn _ Pg (ari\ 

n = — , n N = — , n G = — (461) 

P P P 



satisfy the normalization 

no + n N + n G = 1 . (462) 
Let us define the dimensionless anomalous averages 

a = — , — = a + no ■ (463) 
P P 

Passing to dimensionless quantities, let us use the gas parameter 7, defined in Eq. (365), dimen- 
sionless sound velocity (385), and dimensionless temperature 

« - % ■ < 464 > 
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Then, the equation for the sound velocity (432) reads as 

s 2 = 47r 7 (l - n G + a) . (465) 
For the normal fraction of uncondensed particles, we have 



s 3 ( Q /-oo 1/2 



coth I _ I - , 



cfe^ . (466) 



The glassy fraction is 

■'■■- C + 7s 4/7 (s _ C) 3/7- 



n G = t , z 7777 ttv7 • (467) 



And the superfluid fraction becomes 



4 s 5 f°° (VTT^ - 1) 1 xdx 
n s = 1 n G = — / Y J- . 468 

The Bose-Einstein condensation temperature T c , in the presence of disorder, decreases linearly 
with the increasing disorder strength, as compared to the transition temperature T c °, given by 
Eq. (400), for the system without disorder. The relative transition temperature decrease, for 
C < 1, follows [108] the law 

T — T° It 
ST C = = - — . (469) 

T c ° 9tt v ; 

This is in agreement with Monte Carlo simulations [107]. 

Depending on the relation between the localization length l ioc and the coherence length 

/r|p(r,0)|rfr 
J |p(r,0)|dr 

there can exist three different phases [110]. 

Superfluid phase exists, when the localization length is larger than the coherence length: 

hoc > lcoh ■ (471) 

In this case, the disorder is yet weak and cannot destroy the system coherence. 

Bose glass can occur, when the localization length becomes shorter than the coherence length, 
but yet larger than the mean interparticle distance: 

a < hoc < l coh ■ (472) 

Then a kind of granular condensate can exist, being localized in different spatial regions that are 
separated from each other by the normal nonsuperfluid phase. 

Normal glass appears, when the localization length is shorter that the mean interparticle 
distance: 

hoc < a . (473) 

Therefore no coherence between particles can arise, all of them being localized in separate regions 
of deep random wells. 
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9.5 Bose Glass 



The peculiar phase of the Bose glass is the random mixture of Bose-condensed droplets, localized 
in different spatial regions that are separated from each other by the normal phase. In the 
Bose-condensed regions, the gauge symmetry is locally broken, while in the regions of the normal 
phase, the gauge symmetry is preserved. All these regions are randomly distributed in space 
and it is even possible that they chaotically change their spatial locations. Also, they are not 
necessarily compact and may be ramified having fractal geometry [133,134]. 

Such a randomly mixed system is a particular case of heterophase systems, whose examples 
are ubiquitous in condensed matter physics. In this respect, it is possible to mention paramagnets 
with local magnetic ordering revealing spin waves [135-139], many ferroelectrics [140-143] and 
superconductors [144-149], colossal magnetoresistant materials [150-152], and some other systems 
reviewed in Refs. [54,153-157]. 

The typical features of these heterophase materials are: (i) the embryos of one phase inside 
another are mesoscopic, their characteristic sizes being much larger than the mean interparticle 
distance but shorter than the system length; (ii) the spatial distribution of the embryos, as well 
as their shapes are random; (iii) the system, as a whole, is quasiequilibrium, being either stable, 
or at least metastable, with the lifetime essentially longer than the local equilibration time. 

Such materials, with randomly distributed mesoscopic embryos of one phase inside another 
should be distinguished from systems composed of large stationary domains and from Gibbs mix- 
tures of coexisting macroscopic phases [158] . For the equilibrium macroscopic phases, coexisting 
with each other, one has to consider the interfacial free energy [159]. The notion of interfacial free 
energy arises when one considers uniform macroscopic phases, while the mesoscopic heterophase 
inclusions are nonuniform. For quasiequilibrium mesoscopic embryos, the interfacial regions are 
not well defined, being often ramified and nonequilibrium. Quasiequilibrium embryos of com- 
peting phases are also different from nonequilibrium nuclei arising in kinetic phase transitions 
[160,161]. 

A general approach to treating such random heterophase mixtures has been advanced and 
developed in Refs. [162-171], and summarized in the review articles [54,156,157]. Here, this 
approach is applied for describing the Bose glass. 

Assume that we aim at describing the heterophase mixture of normal uncondensed phase and 
Bose-condensed phase. The condensed phase exists in the form of mesoscopic embryos surrounded 
by the normal uncondensed phase. The effective total volume, occupied by each phase is V v , with 
the index u — 1,2 enumerating the phases. The geometric weight of each phase is 

u>v=y (V = V 1 + V 2 ) , (474) 

where V is the system volume. This weight enjoys the standard probability properties 

W! + w 2 = l, < w v < 1 , (475) 

because of which it can be termed the geometric probability. 

The nonuniform mixture of phases, consisting of mesoscopic embryos, requires the use of 
two types of averages, the statistical averaging over the particle degrees of freedom and the 
configuration averaging over all admissible random spatial phase configurations. Accomplishing 
the configuration averaging over phase configurations [54,156,157] results in the appearance of 
the renormalized Hamiltonian 

H = H 1 ^H 2 , (476) 
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consisting of two terms corresponding to each of the phases. This Hamiltonian is defined on the 

mixture space 

M = U l (g)U 2 . (477) 

By its mathematical structure, this space is the tensor product of the weighted Hilbert spaces 
[54,156,157]. It can be treated as a particular case of a fiber bundle [172]. 
The effective statistical operator for the random mixture becomes 



P = Pi 



P2 , 



(478) 



with p v being the effective statistical operators for each of the phases. 

Thus, after the configuration averaging, we come to the statistical ensemble {p, A4} that 
is the collection of the partial phase ensembles {p^,V. u }. The partial ensembles can be called 
the reduced, or restricted, ensembles, since they are defined on the restricted spaces of micro- 
scopic states, typical of the corresponding phase [173-175]. All details are given in the reviews 
[54,156,157]. 

The Hamiltonians H u are the effective phase Hamiltonians. For the Bose-condensed phase, 



H 1 = H 1 - p N - p l N l - A , 

while for the normal uncondensed phase, 

H 2 = H 2 - p 2 N 2 , 

where N 2 is the number operator for the normal phase. 

The energy Hamiltonian for the Bose-condensed phase reads as 



(479) 



(480) 



= wi J V> f ( 



•0(r) dr + 



+ 



^ t (r)^ t (r')$(r - r')V>(r')V>(r) drdr' , 



(481) 



with the Bogolubov-shifted field operator 

V>(r) =77(r) + ^i(r) . 
For the normal phase, the energy Hamiltonian is 



Ho 



•02 (r) dr + 



+ 



J4( 

j V4(r)V4(r')$(r - r> 2 (r')V 2 (r) drdr' . 



(482) 



The interaction potential here is written in the general form. But, in particular, it can take the 
local form (316). 

In the broken-symmetry phase, we have, as earlier, the densities of condensed and uncon- 
densed particles, respectively, 



Po(r) = Mr) 



Pl (r) = (Vl(r)Vi(r)) . 



(483) 
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And in the normal phase, there is only the density of normal particles 

p 2 (r) = (V4(r)V 2 (r)> ■ (484) 

The numbers of particles in the whole heterophase system are written as follows. The number 
of condensed particles is 

N =w 1 [ p (r) dr . (485) 



Here and below, the integration is over the whole system. The number of uncondensed particles 
is given by the average 

N 1 = (N 1 )=w 1 J Pl (r)dr (486) 

of the number operator 

A\ = w 1 J ^{(r)^i(r) dr . (487) 
The number of normal particles is the average 

N 2 = (N 2 ) = w 2 Jp 2 (r) dr (488) 

of the number operator 

N 2 = w 2 J V4(r)V>2(r) dr . (489) 

The total number of particles in the system is the sum 

N = N + N l + N 2 . (490) 

The related fractions of condensed, n , uncondensed, ni, and normal, n 2 , particles satisfy the 
normalization 

n + ni + n 2 = 1 . (491) 

The system chemical potential is 

u = u n + uini + fx 2 n 2 . (492) 

^From the condition of equilibrium, it follows [12] that 

fi n + fum 

u 2 = ■ = u . (493) 

n + m 

The grand thermodynamic potential is defined as in the previous sections, 

n = -T((\nTre-^ S )) , (494) 

with the double brackets implying the stochastic averaging over the random external potential 
£(r). The geometric weights w u are defined to be the minimizers of the grand potential (494), 
under the normalization condition (475). The latter can be taken into account explicitly by 
introducing the notation 

wi = w , w 2 = 1 — w . (495) 
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Then the minimization of the grand potential (494) implies 

d 2 fl 



dw 



dw 2 



> . 



The first condition gives the equation 



for the weight w, while the second, the stability condition 



(496) 



(497) 



(498) 



Since the right-hand side of inequality (498) is non-negative, the sufficient stability condition is 

,d 2 H, 



dw 2 



> . 



(499) 



Let us use the notation for the single-particle terms of the condensed phase, 

~ V 2 



^t( r ) 



2m 



+ U(r) + £(r) 



ip(r) ) dr 



fio J po(r) dx - pi J pi(r) dr , 



(500) 



and the normal phase, 



$(r) 



M^dr - pi J p 2 (r)rfr, (501) 



respectively. Similarly, we can define the interaction terms for the condensed phase, 

$! = J (ft(r)ft(r')<f>(r - r')V>(r')V>(r))drdr' , (502) 

and the normal phase, 

$ 2 = J (^(r)^(r')<D(r - r> 2 (r / )^ 2 (r))rfrrfr / . (503) 

In the above expressions (500) and (502), it is assumed that the linear in ipi terms are omitted, 
being cancelled by the term A in Hamiltonian (479). 

Then Eq. (497) yields the equation for the geometric weight of the condensed phase 



w 



$2 + K 2 - if 1 

and the stability condition (499) results in the inequality 

$i + $ 2 > . 



(504) 



(505) 



The latter condition shows that the heterophase mixture can exist only for particles with repulsive 
interactions. 
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10 Particle Fluctuations and Stability 



10.1 Stability Conditions 

Fluctuations of observable quantities in statistical systems are characterized by the dispersions 
of self- adjoint operators corresponding to observables. Let A be the operator of an observable 
quantity given by the statistical average (A) of this operator. The fluctuations of this observable 
quantity are quantified by the dispersion 

A 2 (i) = {A 2 } - {A} 2 . (506) 

The fluctuations of an observable quantity are called thermodynamically normal when the 
related dispersion is proportional to N a , with a not larger than one. And the fluctuations are 
termed thermodynamically anomalous if the corresponding dispersion is proportional to N a , with 
a larger than one. In recent literature on Bose systems there has appeared a number of articles 
claiming the occurrence of thermodynamically anomalous fluctuations of the particle number in 
Bose systems everywhere below the transition temperature. 

In the papers [63,176-178] and reviews [5,9,12], it has been explained that the occurrence of 
such anomalous fluctuations contradicts the basic principles of statistical physics and that their 
appearance in some theoretical works is due merely to incorrect calculations. Because of the 
importance of this problem, it is described below, being based on Refs. [5,9,12,63,176-178]. 

The ratio of the operator dispersion to its average value quantifies the intensity of the system 
response to the variation of the considered observable. This response has to be finite in order 
that the system would be stable with respect to the observable-quantity fluctuations. That is, 
this ratio has to satisfy the stability condition [12,63,176-178] 

A 2 ( A) 

< =4^ < oo . (507) 
~ \(A)\ 

This condition must hold for all observables and for any statistical system, including thermody- 
namic limit. For extensive observables, to be considered below, (A) oc N. The limiting ratio 

X (A)=hm ^ (508) 

iV— s>oo iv 

has the meaning of the response function related to the variation of the observable represented 
by the operator A, and can be called fluctuation susceptibility. Therefore, another form of the 
stability condition is 

< x(A) < oo . (509) 

The number of particles is the observable represented by the number operator N. Hence, the 
stability condition with respect to particle fluctuations is 

< x(N) < oo . (510) 

/,From the general relations of statistical mechanics and thermodynamics that can be found 
in almost any course [49,50,59,60,62,79,81,88,89,90], it is easy to show that quantity (508) is 
really proportional to some physical susceptibility. Being interested in particle fluctuations, 
one has to consider the dispersion A 2 (N). The related physical susceptibility is the isothermal 
compressibility. This can be defined in any statistical ensemble, as is shown below. 
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In the canonical ensemble, where the thermodynamic potential is the free energy F = 
F(T, V, N), the isothermal compressibility is given by the derivatives 



/ TN v \ / TN 



In the Gibbs ensemble, with the Gibbs thermodynamic potential G = G(T, P, N), the com 
pressibility is 

1 fd 2 G\ 1 fdV 



^ = -v[^) TN = -v{ap) TN - < 512 > 

And in the grand canonical ensemble, with the grand thermodynamic potential Q = fl(T, V, p), 
the compressibility becomes 

1 fd 2 Q\ 1 fdN\ 

KT --w P W) TV -w P [^) TV - (513) 

Of course, the value of the compressibility does not depend on the used ensemble, provided 
that it is correctly defined as a representative ensemble [54,63,64,71]. The fact that the com- 
pressibility is directly related to particle fluctuations is the most evident in the grand canonical 
ensemble, where 

A 2 (N) , x 

Kt = ^ J . 514 
pTN y ' 

The importance of correctly describing the particle fluctuations is caused by the fact that 
they define not only the compressibility, but also are connected with several other observable 
quantities, such as the hydrodynamic sound velocity St, 

9 1 fdP\ 1 NT , N 

4 = - — = = ^ , 515 

m\dpJ T mpn T m A 2 (N) 

and the central structure factor 

T A 2 (N) 

s (°) = " r '* = ^ = ^ i ' (516 > 

As is seen, the susceptibility x(N) coincides with the structure factor (516). 

Is is worth stressing that all expressions (511) to (516) are exact thermodynamic relations 
that are valid for any stable equilibrium statistical system. 

In stable statistical systems, the compressibility, as well as the structure factor, are finite. 
This is a very well known experimental fact. They can be divergent only at phase transition 
points, where, as is known, the system is unstable. But everywhere outside of transition points, 
all these quantities must be finite. 

10.2 Fluctuation Theorem 

The stability condition (509) is necessary in order that the system would be stable with respect 
to the fluctuations of the observable quantity represented by the operator A. But what can be 
said with regard to an observable represented by a composite operator 
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given by a sum of several self-adjoint operators? How the fluctuations for the total sum of A are 
connected with partial fluctuations for A{! To formulate this question more precisely, let us give 
some definitions. 

Definition: Thermodynamically normal fluctuations 

Fluctuations of an observable quantity, represented by a self-adjoint operator A, are called 
thermodynamically normal if and only if the stability condition (509) holds for this operator. 
Then the related susceptibility (508) is also called thermodynamically normal. 

Definition: Thermodynamically anomalous fluctuations 

Fluctuations of an observable quantity, represented by a self-adjoint operator A, are called 
thermodynamically anomalous if and only if the stability condition (509) does not hold for this 
operator. Then the related susceptibility (508) is also termed thermodynamically anomalous. 

The question of interest is how the total fluctuation susceptibility x(A) is connected with the 
partial fluctuation susceptibilities x(A)? Or, in physical terminology, can it happen that the 
total susceptibility be finite, while some of the partial susceptibilities be infinite? The answer to 
this question is given by the following theorem on fluctuations of composite observables. 

Fluctuation Theorem. (Yukalov [63,177]) 

Let the observable quantity be represented by a composite operator (517) that is a sum of 
self- adjoint operators. Then the dispersion of this operator is 

A2 (e^) = E A2 (^) + E^V A2 (^) A2 A)> ( 518 ) 

where < 1, hence the total fluctuation susceptibility reads as 

* ( E A ) = E XA) + E ^ y/x&MAj) . (519) 

From here it follows that the total fluctuation susceptibility is normal if and only if all partial 
fluctuation susceptibilities are normal. And the total fluctuation susceptibility is anomalous if 
and only if at least one of the partial fluctuation susceptibilities is anomalous. 



10.3 Ideal-Gas Instability 

For illustrative purpose, one often considers the ideal Bose gas. The grand Hamiltonian for 
noninteracting particles is a particular case of Hamiltonian (58), where the energy Hamiltonian 
is 

H = J ^{r) ^- + l?J V>(r) dr . (520) 
Substituting here the Bogolubov shift (48) yields the grand Hamiltonian 

H = J v *(r)(-^- + U- /i ) 77(r) dr + J 4(r) (- ^ + U - pt^j ^i(r) dr . (521) 

The equations of motion become 

d /V 2 \ 

i m T,(r,t)=l-— + U-«>)v(r,t), 
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d / V 2 

t-^(r,t)= l-— + U- f i 1 )^ 1 (r,t). 



(522) 



Let us pass to the uniform gas, when there is no external potential, U — 0. Then, in 
equilibrium, the condensate function is constant, r](r,t) = r] = const. The equation for the 
condensate function gives /xo = 0. 

In the momentum representation, Hamiltonian (521), with U — 0, reduces to 



(523) 



The condition of the condensate existence (15), as well as the Hugenholtz-Pines relation (203), 
result in /ii = for temperatures below the condensation temperature 



2tt 



m 



C(d/2) 



d/2 



(524) 



written here for a <i-dimensional space. Expression (524) shows that positive T c does not exist 
for d — 1, since C(l/2) = —1.460 and that T c = for d = 2, since £(1) = oo. Positive T c exists 
only for d > 2. 

For the number operator iV = N + Ni, taking into account that cov(iV ,iV 1 ) = and 
A 2 (iV ) = 0, one finds 

A 2 (iV) = A 2 (iV 1 ) . (525) 

Let us emphasize that the condensate fraction does not fluctuate at all and that the total 
fluctuations are caused solely by the uncondensed particles. The number operator of the latter 
is 

N, = f Vl(r)Vi(r) dv = J2 a W ■ (526) 

Invoking the commutation relations and the Wick theorem, one has 

(a\a k ala p ) = (4aja fc a p ) + 5 kp n k , 
(a\ala k a p ) = n k n p + 5 kp n 2 k , 



where the momentum distribution is 

n k = (a\a k ) = 

This leads to 

^i 2 : 



exp 



2m 



-\ -i 



^1+5^^(1 + ^) • 
Therefore particle fluctuations are characterized by the dispersion 



(527) 



Remark. In some works, the authors forget that Bose-Einstein condensation necessarily 
requires broken gauge symmetry. Forgetting this, one extends the sum in Eq. (526) to k — 0. 
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Then, separating the term with k — 0, one gets the condensate fluctuations described by the 
dispersion A 2 (iVo) proportional to Nq. One blames the grand canonical ensemble to be guilty for 
this unreasonable result, naming this "grand canonical catastrophe". However, as is clear, there 
is no any catastrophe here and not the grand ensemble is guilty, but the authors doing incorrect 
calculations. One should not forget that, if the gauge symmetry is not broken, then N = 0. 

Summing the right-hand side of Eq. (527) yields 

A 2 (iV 1 )= (^) 2 ^ 4/3 - (528) 
This gives the fluctuation susceptibility 

!!^_] N 1 ' 3 = oo. (529) 

7T J 

Consequently, the stability condition (510) does not hold. This means that the ideal uniform 
Bose-condensed gas is not stable. It is a pathological object that cannot exist in reality. 



10.4 Trapped Atoms 

But maybe the ideal Bose-condensed gas could be stabilized by confining it inside a trap formed 
by a trapping external potential? A general expression for such a trapping potential is given by 
the power-law form 



2 

a=l 



a 

In 



(530) 



which is written here in the <i-dimensional space. The trapping frequency u a and the trapping 
length l a are connected by the relations 

1,1 . , 

uj a = — - , l a = = . (531) 

ml* ^mu a 

It is also useful to introduce the effective frequency and effective length by the geometric averages 

]>«) =4- io "(n/"j (532) 

Let us define the confining dimension [51] 

a=l 

Passing from the trapping potential to the uniform case implies the limits 

L d 

n a ->■ oo , l ->■ - , J J 2l a ->■ L d , 



2 

a=l 



where L is the linear size of the system volume V = L d . As a result, s — >■ d/2, that is, s becomes 
semi-dimension. 
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Bose-Einstein condensation of the ideal Bose gas in the trapping potential (530) can be 
described employing the generalized quasi-classical approximation [51]. The condensation tem- 
perature reads as 



bN 



where we use the notation 



6 = ^rllf7T 



1/s 



l/2+l/n« 



(534) 



JUL r(l + l/n a ) 



and introduce the generalized Bose function 



9s(z) = 



1 



00 zu"- 1 



du , 



in which the integration is limited from below by the value 

Uo =2T- 



(535) 



(536) 



Recall that in the standard Bose function, the integration starts from zero. 

The value <? s (l) of the generalized function (535) is finite for all s on the complex plane, 
since T(s) ^ 0, so that 1/T(s) is an entire function. But T(s) can be negative for s < 0, e.g., 
it is negative in the interval — 1 < s < 0. Therefore, <7 S (1) is positive and finite for all s > 0. 
Contrary to this, the standard Bose function would diverge for s < 1, and there would be no 
finite condensation temperatures for these s. While, in the case of the generalized function (535), 
finite condensation temperatures formally exist for any positive s. Below T c , and for s > 0, the 
condensate fraction is 

n = 1 - (T < T c ) . (537) 

2 and 



The most often studied trapping potential is the harmonic potential, for which n a 
s = d. Then the condensation temperatures are 



Nu 



Tc ln(2JV) 



UJ 



N 



1/d 



(d=l), 



(d > 2) . 



(538) 



The condensation temperature (534) is finite for any finite N. But it is necessary to check 
whether it is finite in thermodynamic limit, when iV — > oo. For confined systems, the effective 
thermodynamic limit is defined [51] in Eq. (13). As an extensive observable, we can take the 
internal energy that, in the present case, below T c , is 



E N = - g 1+s (l)T 1+s . 



Then, the thermodynamic limit (13) reads as 

N — > oo , E N — > oo , 



En 
N 



— > const 



(539) 



(540) 
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The value g 1+s is finite for iV — > oo at all s > 0. Hence, Eq. (540) can be rewritten as the limit 

iV->oo, b^O , bN -> const . (541) 

For the equipower traps, for which n a = n, the effective thermodynamic limit (541) takes the 
form 

N ->■ oo , cu ->■ , A^Wg -> const . (542) 

Considering the thermodynamic limit for the condensation temperature (534), we have to 
take into account that the generalized function (535) yields 

1 /2T\ 1_S 
(1 -s)T(s) \u J 

g s (l) = In (s = 1) . (543) 

Consequently, for the condensation temperature (534), as iV — > oo, we find 

T < a JVO^IA ^ «>< " < 1) , 

T ^kT^° = 

T c ->■ const (s > 1) . (544) 

Therefore, finite condensation temperatures exist only for s > 1. This implies that for harmonic 
traps, for which s = d, the finite condensation temperature occurs only for d > 2. Bose-Einstein 
condensation cannot happen in one-dimensional harmonic traps at finite temperature. 

But this is not yet the whole story. As we know from Sec. 10.3, a finite condensation 
temperature can formally occur, however, the condensed system in reality is unstable, thus, 
cannot exist. To check the stability, it is necessary to consider the system susceptibilities. Specific 
heat for the Bose-condensed trapped gas is finite at all temperatures, displaying a jump at the 
transition point [51]. We need to consider the isothermic compressibility (514) that shows the 
system response with respect to particle fluctuations. The dispersion for the number operators 
behaves as 

A 2 (iV ) = , A 2 (iV) = A 2 (7Vi) . (545) 
It is convenient to introduce the finite- N susceptibility 

A 2 (iV) , x 

Xn = — ^ , (546) 



whose limit 



lim xn = x( N ) 

N— >oo 



yields the susceptibility defined in Eq. (508). Below T c , we obtain 



& = ' (547) 



75 



Susceptibility (547) is negative for s < 1 and does not satisfy the stability condition (510). 
For the values of s > 1, we have 

2 ( T x " 

Xn 



N \ T t / 

\v \77~TF7 7t(—) (L) ■ (K-s<2). 



(2 - s)c(s)r(s - 1) V^o / V^e 



1 /T\ 2 , /2T\ 



«"=i?(3' (s>2 »- < 548 » 

For asymptotically large JV, we get 

Xn oc N (s = 1) , 

Xv oc jV( 2 " s )/ s (l< s < 2) , 
XvoclniV (s = 2), 

Xv oc const (s > 2) . (549) 

This shows that the trapped Bose gas is stable only for s > 2, when the stability condition (510) 
is satisfied, that is, when 

d / 

s 

a=l 

In particular, for harmonic traps, for which s — d and b = u d , one finds 



^U±L > 2 . (550) 



\UJq J \ UJq J 

6C(3) Uc 



7T 2 (T 
XN = TTTTTTT I TfT ) (d = 3) . 



For large N, this gives 



Xv oc N (d=l). 

XnocIiiN (d = 2). 

Xn oc const (d — 3) . 

Thus, the Bose-condensed gas in a harmonic trap is stable only in the three-dimensional space, 
d = 3. 

The above analysis demonstrates that confining the ideal Bose gas in a trap may stabilize it, 
which, however, depends on the confining dimension s, defined in Eq. (533). The occurrence of a 
formal expression for the critical temperature T c is not yet sufficient for claiming the possibility of 
Bose-Einstein condensation in a trapped gas, but it is also necessary to check the system stability. 
For example, in the case of the power-law trapping potentials, the condensation temperature 
formally exists for s > 1. But the trapped condensed gas can be stable only for s > 2. One- 
and two-dimensional harmonic traps are not able to stabilize the condensate. Only the three- 
dimensional harmonic trap is able to host the ideal Bose-condensed gas. 
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10.5 Interacting Systems 



Ideal gases are, actually, rather artificial objects, since there always exist particle interactions, 
though, maybe, weak. Now we pass to studying the stability of interacting systems. 

To a great surprise, there have been published many papers, in which the authors claim that 
interacting Bose-condensed systems, both uniform as well as trapped, exhibit thermodynamically 
anomalous particle fluctuations of the same kind as the ideal Bose gas, with the number-operator 
dispersion (528). By the Bogolubov theorem, one always has A 2 (iV ) = 0, hence, A 2 (N) = 
A 2 (A r 1 ). Then the thermodynamically anomalous dispersion A 2 (iV) oc iV 4 / 3 would lead to xn 
iV 1 / 3 and to the divergence of x(N) — > oo. In that case, the stability condition (510) is not 
satisfied, and the behavior of all physical quantities would be rather wild. Then the isothermal 
compressibility (514) would diverge, the sound velocity (515) would be zero, and the structure 
factor (516) would be infinite. That is, the system would be absolutely unstable. 

Moreover, this would mean that any system with spontaneously broken gauge symmetry 
would not exist. Clearly, such a strange conclusion would contradict all known experiments 
observing Bose-Einstein condensed trapped gases. Superfluid helium is also the system with 
broken gauge symmetry, hence, it also would not be able to exist, which is evidently absurd. 

In Refs. [5,9,12,63,176-178], it has been explained that the occurrence, in some works, of 
thermodynamically anomalous fluctuations is caused by incorrect calculations. One calculates the 
dispersion A 2 (iV\) invoking the Bogolubov approximation that is a second-order approximation 
with respect to the operators of uncondensed particles. But the expression Nf is of fourth 
order with respect to these operators. Calculating the fourth-order terms in the second-order 
approximation, strictly speaking, is not self-consistent and can lead to unreasonable results, such 
as the occurrence of thermodynamically anomalous particle fluctuations. 

The correct calculation of the dispersion A 2 (iV) and, hence, of susceptibility (546), can be 
done as follows. From the definition of the particle dispersion A 2 (iV), one has the exact expression 



which is valid for any system, whether uniform or nonuniform, equilibrium or not [176-178]. Here, 



is the pair correlation function. 

In the HFB approximation, analogously to the Bogolubov approximation, one has to retain 
in Eq. (551) the terms up to the second-order with respect to the operators of uncondensed 
particles. For nonuniform systems, one can employ the local-density approximation of Sec. 7.4. 
Then Eq. (551) reduces to 




(551) 



p(r) = po(r) +pi(r) 



is the total particle density and 



(^t(r)^t(rO^(r')V>(r)) 
p(r)p(r') 




(552) 



Using the formulas of Sec. 7.4 gives 




(553) 
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The same Eq. (553) represents the structure factor (516). Equation (515), defining the hydro- 
dynamic sound velocity, leads to 

-i 

(554) 

And the isothermal compressibility (514) becomes 

provided the average density p is defined. 

For a uniform system, the above formulas reduce to 

A 2 (N) T 1 

Xn = — ^— = S(0) = — -, s T = c, k t = -. (556) 

JM mcr mpc z 

It is important to emphasize the necessity of taking into account the gauge symmetry breaking 
in the above calculations. If the symmetry would not be broken, or if the anomalous average 
a would be omitted, one would get the divergence of expressions (553) and (555), which would 
mean the system instability [179]. 

In some works on particle fluctuations, one also makes the following mistake. One writes 
that, in the canonical ensemble, the condensate fluctuations are given by A(A r ) that is equal 
to A(A r 1 ), and one calculates the latter in the second quantization representation. However, 
this representation uses the field operators defined on the Fock space and, by construction, it 
is introduced for the grand canonical ensemble. So, A(A r i) has nothing to do with condensate 
fluctuations that, by the Bogolubov theorem correspond to A(A r ) = 0. 

In this way, correct calculations lead to no anomalous thermodynamic particle fluctuations. 
The latter arise only in incorrect calculations. There are no anomalous fluctuations neither in 
correctly employed Bogolubov or HFB approximations [5,9,12,63,176-179] nor in the renormal- 
ization group approach [180]. 

If thermodynamically anomalous fluctuations would not be caused by calculational defects, 
but would be real, then not merely equilibrium Bose-condensed gas and superfluid helium would 
not exist, but the situation would be even more dramatic. This is because the systems with gauge 
symmetry U(l) are just a particular case of systems with continuous symmetry, all such systems 
having general properties connected with their continuous symmetry and the symmetry breaking 
[181]. Therefore all such systems exhibiting thermodynamically anomalous fluctuations would 
not exist. We mean here only equilibrium statistical systems, since nonequilibrium systems can 
possess strong fluctuations making them unstable [182-184]. 

For example, many magnetic systems exhibit continuous symmetry connected with spin ro- 
tation. The appearance of magnetic order in such magnetic systems implies the spontaneous 
breaking of the spin rotational symmetry. If the continuous symmetry breaking would lead to 
the appearance of thermodynamically anomalous fluctuations of the order parameter, then, in 
magnetic systems, this would mean the occurrence of thermodynamically anomalous magnetic 
susceptibility, hence, instability. Then there would be no stable equilibrium magnetic systems 
with continuous symmetry breaking, which is again absurd. 

To show that the spontaneous breaking of the spin-rotation symmetry does not lead to ther- 
modynamically anomalous magnetic fluctuations [12], let us consider the Hesenberg model, with 
the Hamiltonian 



1 

N 



P(r) 
c 2 (r) 



dr 
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in which Sj is a spin operator on the i-lattice site, Jij is an exchange interaction potential, and B 
is an external magnetic field. The Hamiltonian enjoys the spin rotation symmetry in the absence 
of the external field B. 

The Gibbs potential is defined as 

G — -Tin Tre~ pA = G(T, N, B) . (558) 

The system magnetic moment 

Ms m = - { m )s &> (559) 

can be represented as the average of the magnetic-moment operator 

M = " M = 5>- ( 56 °) 

i 

The magnetic susceptibility tensor is given by the elements 

1 dM p _ 1 d 2 G 
XaP ~ N dB a ~ N OB a dB p ■ [ } 

Direct calculations yield 

Xa p = ^cav{M a ,Mp), (562) 
which shows that the diagonal elements 

Xaa = NT (563) 

are expressed through the dispersion of the components of the magnetic-moment operator (560). 

The HFB approximation for the field operators of Bose systems is equivalent to the mean-field 
approximation for the spin operators of magnetic systems. Therefore, it is reasonable to resort 
here to the mean-field approximation, although the results are qualitatively the same if we invoke 
more elaborate techniques. In the mean-field approximation, Hamiltonian (557) reads as 

H = - H • S i + N J($i) 2 , (564) 

i 

in which the notation is used for the effective field 

H = 2 J(Si) + B (565) 

and the effective interaction 

J 

X 



J = jj E J v ■ ( 566 ) 

For concreteness, let us consider spin one-half. Then the Gibbs potential (558) becomes 



G = -NT In 



2 cosh 



Ho 
2T 



+ NJ(S l )\ (567) 
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where the ideality of the lattice is implied and 



Ho = \n\ = JY, H2 «- 



The average spin is defined by the extremization condition 

dG 



0, 



which is equivalent to the equation 

{ l) ~ N dB ■ 

As a result, one finds 

< S -> = 2^ ^ (I 
This defines the magnetic susceptibility (561) as 

X aP = ^ (Sf) • 

Let us define the order parameter 

V = 2\(Si)\. 

In view of Eq. (570), this reads as 

H 



ii = tanh I — - 
I 2T 



Susceptibility (571) takes the form 



V 

: 11 ■ I,, ., ,\-. ,/ \ ( 1 >r I 



Directing the external magnetic field along the axis z, so that 

B x = B y = , B z = h , 

yields 

H x = H y = , H z = 2 J (Sf) + 5 2 , 

which can be rewritten as 

-Ha — SazHo , H = H z = Jl] + h . 

The average spin components become 

(5?) = (Sf) = , (Sf) = I . 
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The susceptibility tensor (574) leads to the equation 



Xa/3 = ^(<W + 2JXc 



?7 ( 1 - rf 

— + °0z 



(578) 



rj = tanh ( ^±^1 . (579) 



_H "* \ 2T H 

with the order parameter 

, (J 

n = tanh I 

22 

Equation (578) shows that the nondiagonal elements are zero: 

Xxy = Xxz = Xyz = , (580) 
while the diagonal elements give the transverse components 



and the longitudinal component 



Xxx = Xyy = 7^ ( 581 ) 



1 -rf 

Xzz = 2[2T- J(l-v 2 )] ' (582) 



The latter, with the notation for the critical temperature T c = J/2, can be represented as 

l-rf 

Xzz = 4[T-T c (l-7? 2 )] • (583) 
At low temperature, and h — > 0, the order parameter (579) behaves as 

V - 1 - 2 exp (- ^ (T « T c ) (584) 

and susceptibility (583), as 

Xzz * \ exp (- ^ (T « T c ) . (585) 
At high temperature, and h — > 0, susceptibility (583) acquires the Curie- Weiss law 

' (T > T c ) . (586) 



4(T - T c 



The latter susceptibility diverges at the critical point T c . However, this divergence has nothing 
to do with the thermodynamically anomalous behavior, since this is the divergence with respect 
to temperature T, but not with respect to the number of particles N. In addition, the phase 
transition point is the point of system instability, where the system becomes nonequilibrium and 
fluctuations have right to infinitely rise. 

One introduces the transverse susceptibility 

X± = Xxx = —^f 1 (587) 
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and the longitudinal susceptibility 

X|| = x« = , (588) 

where relation (563) is taken into account. In view of Eqs. (581) and (582), one finds the 
dispersions for the magnetic- moment operator (560), characterizing the transverse fluctuations, 

A 2 (M X ) riT , . 

Tx± = ir, (589) 



N ^ 2h 



and the longitudinal fluctuations, 



^v" = TX|1 = 4[T-T C (1-^)] ' (590) 

of the system magnetic moment. 

The longitudinal fluctuations are always thermodynamically normal, if calculated in a self- 
consistent way, as it should be, according to the stability condition (509). In some papers, one 
finds thermodynamically anomalous magnetic fluctuations of the same type as for Bose systems, 
with A(M Z )/N oc iV 1 / 3 . But, as has b een explained in Ref. [5], this is due to the same mistake 
as one does when dealing with Bose systems. One approximates Hamiltonian (557) by a second- 
order form, with respect to small deviations from the average magnetic moment. And then, one 
considers the fourth-order form calculating the dispersion of M z . Going outside of the region of 
applicability of the chosen approximation leads to the appearance of meaningless results. But 
self-consistent calculations, as is shown above, always give normal longitudinal fluctuations. 

The transverse fluctuations are known [185] to be much larger than the longitudinal ones. 
Formally, Eq. (589) diverges when h — > 0. This, however, does not make the transverse magnetic 
fluctuations thermodynamically anomalous. To be thermodynamically anomalous, expression 
(589) should diverge with respect to the number of particles N, or, what is the same, with 
respect to the system volume V. But here, it is the divergence with respect to h. 

Moreover, one should not forget that below the transition temperature T c the spin-rotation 
symmetry is broken. The symmetry breaking is described by switching on a small external mag- 
netic field h 7^ 0. But then Eq. (589) is finite. Switching off this field restores the symmetry, as a 
result of which Eq. (589) would diverge, similarly to how the compressibility of a Bose-condensed 
system would diverge being incorrectly calculated without the gauge symmetry breaking. There- 
fore, as soon as the spin-rotation symmetry has been broken, when h ^ 0, all fluctuations are 
thermodynamically normal. And above T c , where the symmetry is not broken, one has rj = h/2T, 
hence A 2 (M X )/N = 1/4, which is again finite for any h. 

When some symmetry in a system is broken, the mathematically correct definition of statis- 
tical averages is understood in the sense of the Bogolubov quasiaverages [16]. Then, as is well 
known, one has, first, to accomplish the thermodynamic limit, with N — > oo and, only after 
this, to consider the limit h — > 0. In that sense, there is no any thermodynamically anomalous 
fluctuations. 

Note that real magnetic systems always possess magnetic anisotropy. This can be small, 
but never exactly zero, which corresponds to the presence of a finite h. Consequently, in real 
equilibrium magnetic systems, there are no thermodynamically anomalous fluctuations. And 
there are no thermodynamically anomalous fluctuations in any equilibrium system with the 
spontaneous breaking of any continuous symmetry. In the other case, such a system would be 
unstable and could not be in equilibrium. 
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11 Nonground- State Condensates 



11.1 Coherent Modes 



First of all, it is necessary to concretize what is meant under nonground-state condensates. 
The stationary equation (335) for the condensate function can be treated as an eigenproblem. 
Generally, an eigenproblem can yield a spectrum of possible eigenvalues and a set of the related 
eigenfunctions. So, generally, the eigenproblem, corresponding to Eq. (335), can be represented 
in the form 

V 2 

+ U(r) 



2m 



Vn(r) + 



+ $ [|77 n (r)|V(r) + 2pi(r)77 n (r)+<7i(r)7 7 ;(r)+e(r)] = E n r] n (r) , 
in which the minimal eigenvalue defines the chemical potential 

Ho = min E n . 



(591) 



(592) 



When E n = /xq, Eq. (591) corresponds to the standard ground-state Bose-Einstein condensate, 
while, for higher eigenvalues E n , this equation corresponds to nonground-state condensates. The 
values of pi(r), o"i(r), and £(r) depend on the index n, but for short, this dependence is not 
shown explicitly. 

The condensate function describes the coherent part of the system. In the limit of asymptot- 
ically weak interactions ($o — > 0) and low temperature (T — > 0), when the whole system is in 
the coherent state, Eq. (591) reduces to the nonlinear Scrodinger equation 



V 2 



r) n (r) + $ |r/ n (r)| 2 ?7„(r) = E n rj n (r) . 



(593) 



In the particular case, for E n = /i , it is called the Gross-Pitaevskii equation. 

The condensate function r) n (r) is normalized to the number of condensed particles, as in Eq. 
(52). It is convenient to introduce the function y? n (r) by the relation 



N (p n {r) , 



so that (p n (r) be normalized to one, 



Then Eq. (591) transforms into 



|v?„(r)| 2 dr = 1 . 



(594) 



(595) 



2m 



+ U(r) 



No|<Pn(r)lVn(r) + 2pi(r)<p n (r) + <7i(r)<p*(r) 



E n ip n (r) , 



while Eq. (593), into 



2m 



+ U(r) 



(f n (r) + $o-^o|Vn(r)|Vn(r) = E n ^ n (r) 



(596) 



(597) 
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The solutions to Eqs. (591) and (596) define the coherent modes in the general case [66] and 
Eqs. (593) and (597), the coherent modes for asymptotically weak interactions and temperature 
[186]. These coherent modes, first introduced in Ref. [186], correspond to nonground-state con- 
densates. The properties of such coherent modes and the methods of their generation have been 
studied in a series of papers [66, 186-215]. A dipole coherent mode was excited in experiment 
[216]. These coherent modes are also called topological, since the nonground-state condensates, 
corresponding to different coherent modes, describe particle densities with different spatial topol- 
ogy 



11.2 Trap Modulation 

There are several requirements that are necessary for creating a nonground-state condensate. 
First of all, it is clear that such a condensate cannot be equilibrium. Hence, its creation requires 
the action of external time- dependent fields. Second, the system of Bose particles has to possess 
a discrete spectrum in order that it would be possible to distinguish the usual ground-state Bose- 
Einstein condensate from a nonground-state condensate. This means that the system is to be 
placed inside a trapping potential. And, third, to transfer particles from their ground-state to 
a chosen excited state, it is necessary, either to employ a resonant field or to use rather strong 
pumping. 

A straightforward way of imposing external alternating fields is by modulating the trapping 
potential. Let the confining potential be composed of two parts, 

U(r,t) = U(r) + V(r,t) , (598) 

in which the first term is a trapping potential and the second term 

V(r,t) = V 1 (r)cosout + V 2 (r)smout (599) 

realizes the modulation of this potential with frequency ou. 

There exists one limitation on the spatial dependence of the modulated trapping potential 
(598). In Refs. [206,207], the shape- conservation theorem has been proved, showing that the trap 
modulation moves the whole condensate without changing its shape if and only if the trapping 
potential U (r) is harmonic, while the modulation term (599) is linear with respect to the spatial 
variables. In that case, the trap modulation would not be able to produce excited coherent 
modes. So, to generate these modes, one has to avoid this particular case of spatial dependence. 

Suppose that at the initial time t — the system has been completely condensed, being in 
the energy state E — /i . Then, to transfer the system to an energy state E n , one has to use the 
alternating field with a frequency ou close to the transition frequency u n = E n — /i . Under this 
resonance condition 

— < 1 (Alu = ou - ou n ) , (600) 
ou 

it is sufficient to invoke the pumping fields of small amplitudes. 

The time-dependent equation (323) for the condensate function, after the substitution of the 
relation 

r}(r,t) = VNotp(T,t), (601) 

similar to Eq. (594), transforms into 

V 2 1 

+ U(v,t)-fi (p{r,t) + 



^V?M) 



2m 
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+ $0 



iV p (r, t)<p{r, t) + 2pi(r, t)tp(r, t) + <n(r, *)^*(r, t) + 



(602) 



We can look for the solution to this equation represented [66,186,187,198] as an expansion 
over the coherent modes, 

V9(r,t) = ^C n (t)v9 n (r)e-^, (603) 

n 

so that the coefficient function C n (t) be slow as compared to the fast oscillating exponential 
functions: 

(604) 



1 



dt 



< 1 . 



Let us introduce the matrix elements corresponding to particle interactions, 



a mn = %N J |y? m (r)| 2 [2|v?„(r)| 2 - |y? m (r)| 2 ] dr , 



and to the action of the modulating field, 



V4(r) [V^-iV^ipni^dr. 



(605) 



(606) 



Also, let us define the expression 



e n (t) = a r 



% J <(r) 



2p < r\v)<p n {v) - 2p 1 (v,t) iPn (v) + ,r(r)^(r) + 



(«), 



dr , 



(607) 



in which the functions with the upper index n correspond to the stationary solutions characterized 
by the condensate function y?n( r ) and where 



N J \<p n ( 



r)\ 4 dr 



The trap modulation produces not only the required coherent mode but it also destroys 
the condensate by transferring particles from the condensate to the fraction of uncondensed 
particles. Therefore, the generation of the coherent mode can be effectively done only during a 
finite depletion time t dep , when the transfer from the condensate to the uncondensed fraction is 
yet negligible. During this time, the variation of quantity (607) is small, such that 



t de n 
e n dt 



< 1 (t < t dep ) . 



It is convenient to make the change 

C n (t) = c n (t) exp[-ie n (i)t] , 



(608) 



(609) 



in which e n = e n (t) is treated as a slow function of time, in the sense of inequality (608). We 
may notice that in the limit of a completely coherent system, when the fraction of uncondensed 
particles is negligibly small, Eq. (607) does not depend on time. 

Then, we substitute expansion (603) into Eq. (602), employ the above notations, and invoke 
the averaging techniques [114,117,119,121,217-219], based on the existence of different time scales 
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[220,221]. As initial conditions, we assume nonzero c n (0) and c (0), while all other coefficient 
functions c,(0) = for j ^ 0, n. This procedure yields [66,186,195] the equations 

. (ICq I 1 2 1 q iAujt 

1 —J- = ttOn|C n | C + - P0nC n e , 

at 2 

1 IF = a «°l c °| 2c « + \ ^onCoe-^* • (610) 
Solving these equations gives the fractional mode populations 

p n (t) = \ Cn (t)\ 2 . (611) 

As a concluding remark to this section, it is worth emphasizing that expansion (603) corre- 
sponds to the diabatic representation [66,215], and one should not confuse it with the adiabatic 
representation [222], which is not suitable for the studied resonant process. 



11.3 Interaction Modulation 

Another way of exciting the cloud of particles confined inside a trap is by varying the particle 
interactions by means of the Feshbach-resonance techniques [1-3,5,34,223]. This method can also 
be used for generating the coherent modes, as has been mentioned in Refs. [206,207,210], and 
analyzed in detail in Refs. [214,215]. 

Let the scattering length be modulated so that the particle interaction becomes time-dependent 
according to the law 

$(t) = $ + cos(wt) + $ 2 sm(ut) . (612) 
Following the same procedure as in the case of the trap modulation and introducing the notation 

lm = iVo($i - i$ 2 ) f </?o(r)l</?m(r)|Vn(r) dv , (613) 

in which n is fixed and m — 0, n, we get the equations 

• ^ C I 1 2 i / I |2 , 1 I |2 1 iAujt i 1 * * 2 -iAuit 

1 ~d~t = a ° n °° w 1 c ° I +2 7ri ' Cri l J c " e +2 7 o c « c o e 

« = "noM c n + I 7 n |c„| +-7ol c o| J c e + - 7nC c n e . (614) 

Both these ways of modulating either the trapping potential or particle interactions can be 
used for generating excited coherent modes. 

Nonequilibrium systems with the generated coherent modes, representing nonground-state 
condensates, possess a variety of interesting properties. We can mention the following effects: 
interference patterns and interference currents [194,195,198], mode locking [186,198,199], dynam- 
ical phase transitions and critical phenomena [190,194,195,198], chaotic motion [206,207], atomic 
squeezing [198,201,202], Ramsey fringes [211-213], and entanglement production [224-226] that 
can be quantified by a general measure of entanglement production [227-229]. 

The above-mentioned effects can be realized by resonant alternating fields of rather low am- 
plitudes. When increasing the amplitude of the pumping field, it becomes feasible to generate 
the excited coherent modes with the frequencies of the alternating fields, which are not exactly 
in resonance with the transition frequencies. Thus, the transition between the coherent modes, 
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characterized by the transition frequency Ui 2 , can be done by means of the harmonic generation 
and parametric conversion [206,207]. 

Harmonic generation occurs, when the driving frequency u satisfies the condition 

nu = uj\2 (n = 1, 2, . . .) . (615) 

Parametric conversion requires the use of two alternating fields, with the driving frequencies oj\ 
and u>2, such that 

u>i ± u 2 = ■ (616) 
In the case of two pumping fields, there exists the combined resonance under the condition 

n l co l + n 2 w 2 = oj 12 (ni = ±1, ±2, . . .) . (617) 

And, generally, the application of several external alternating fields, with the driving frequencies 
u>i, can generate coherent modes under the condition of generalized resonance, when 

^UiUi = ui2 (rii = ±1,±2, . . .) . (618) 

i 

The amplitudes of external alternating fields can be made arbitrarily strong. Therefore, all 
these effects can be realized in experiments. The particle interactions can also be varied in a 
wide range. For instance, employing the Feshbach resonance techniques, it is possible to tune 
the interactions of 7 Li atoms over seven orders of magnitude [230]. Hence, the generation of 
nonground-state condensates can be done by the interaction modulation as well. 

11.4 Turbulent Superfluid 

As follows from the previous sections, increasing the modulation amplitude results in the genera- 
tion of more and more coherent modes, whose excitation becomes more and more easy, especially 
when several alternating fields are involved. This is because it is sufficient that the frequencies 
of the modulating fields be such that one of the above conditions be approximately satisfied. As 
soon as this happens, the related coherent modes become excited. Intensive field modulation 
generates simultaneously several coherent modes. 

When alternating fields are applied creating an oscillating anisotropy, with local rotation 
moments, then the prevailing coherent modes will be quantum vortices. An important feature of 
the vortex creation by means of the anisotropic trap modulation, contrary to the vortex creation 
by means of rotation, is the generation of vortices as well as antivortices, that is, the generation 
of the vortices with opposite rotation velocities. The oppositely rotating vortices repel each 
other and diffuse in space, separating from each other. At the beginning, when the amplitude 
modulation is not yet too strong, there should arise just a small number of vortices having the 
standard properties [1-3,231,232], except that vortices and antivortices both are present. In the 
case, when the whole system is uniformly rotated, the increased rotation frequency induces a 
vortex lattice [231,232]. Contrary to this, when the trapped system is subject to the action 
of alternating fields, nonuniformly and anisotropically shaking the trapped particle cloud, the 
created vortices possess different axes of rotation and different rotation velocities. Therefore 
no vortex lattice is possible. Then, increasing the amplitude of the alternating fields produces 
a large number of vortices with randomly distributed vorticities. Such a tangle of quantized 
vortices forms what is called quantum turbulence, and the whole system is said to be in the state 
of turbulent superfluid. 
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The problem of turbulent superfluid has been addressed in a number of works. The related 
literature has been reviewed in articles [233-237] (see also recent Refs. [238-240]). There are 
plenty of experiments observing quantum turbulence in liquid 3 He and 4 He. Quantum turbu- 
lence in trapped gases has also been observed [241]. The description of turbulent superfluids as 
continuous vortex mixtures has been advanced [237]. 



11.5 Heterophase Fluid 

Increasing further the amplitude of the alternating field breaks the turbulent superfluid into 
spatially separated pieces, with Bose-condensed droplets separated by normal, nonsuperfluid, 
spatial regions. Such a state reminds the Bose glass, or granular condensate, considered in Sec. 
9.5. But now it is a highly nonequilibrium state. This state is analogous to heterophase mixtures 
consisting of several randomly intermixed phases [54] . Thence, it is called heterophase fluid. 

An external modulating field acts on the system similarly to the action of a spatial random 
potential [110], such as treated in Sec. 9. The possibility of mapping the system with a time- 
dependent modulation to the system with a spatially random potential is very important, since it 
allows us to understand the behavior of modulated nonequilibrium systems by comparing them 
with equilibrium random systems. The proof of this mapping is as follows. 

Let the system Hamiltonian 

H{t) =H + V(t) (619) 
consist of the usual term H , containing no time-dependent fields, and a term 

V(t) = J ft(r)V{r,t)ip(r) dr , (620) 

with an external potential depending on time. The characteristic variation time t mo d of the 
modulating potential V(r,t) is assumed to be much longer than the local-equilibrium time ti oc , 
but much shorter than the time of experiment t exp , 

hoc <^ tmod texp ■ (621) 



1 ptexp 

= N 



dt . (622) 



The modulating potential pumps energy into the system that can be associated with the 
effective temperature 

If the pumping potential is periodically oscillating with a frequency oj and period t, mo d = 2ty/uj, 
then 

9V(t) ~, . 2ir ~ . 
= t»V(t) = — V(t) . 

''mod 

In this case, the effective temperature (622) is 

T* = -f— / | (V(t)) | dt . (623) 



Nt mo< i Jq 

Denoting the amplitude of the modulating potential V(r,t) as V mo d, we have 

| (V(t)) | » NV mod . 
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Therefore, the effective temperature (623) becomes 

T* 2tt h^IL Vmod . (624) 

tmod 

One may notice that the effective temperature depends on t exp , though this dependence is week, 
in the sense that 



tmod 



dT* 



dt 



Under the slow modulation, such that t mo ^ ^> ti oc , the system, at each moment of time, is in 
quasiequilibrium. Consequently, one can define the local in time thermodynamic potential 

n(t) = -T* \nTrexp{-/3*H(t)} , (625) 

where f3* = 1/T*. Because t exp 3> t mod , we are interested not in the local potential (625) but in 
the coarse-grained potential 

j Pernod 

D, = / n(t) dt , (626) 

tmod Jo 

averaged over oscillations that are fast as compared to t exp . 

At each moment of time t the potential V(r,t) describes a spatial potential. This can be 
characterized by the relation 

V(r,f)=£(r), (627) 

which defines the functional 

t = t[t(r)] . (628) 

Equations (627) and (628) symbolize the fact that for each time t there corresponds a spatial 
potential £(r) and vice versa, a potential £(r) is ascribed to time t. The relation between the 
interval [0,t mO( i] and the topological space {£( r )}> without much loss of generality, can be taken 
as homeomorphic. 

The variation of time is equivalent to the variation of the spatial potential, so that 

stm\ 

With relation (628), Hamiltonian (619) becomes the functional 

H[£(t)] = H + J ^{v)i{v)^{v) dr (629) 

of the spatial field £(r). Therefore, the averaged thermodynamic potential (626) takes the form 

Q = -T* J lnTrexp{-/r#[£(r)]} Df(r). (630) 

The latter is equivalent to the thermodynamic potential of an equilibrium system in a random 
external field. 

If the external alternating field has an amplitude V mod and the whole trapped system is 
subject to the modulation, then the modulation amplitude V mod plays the role of the correlation 
amplitude V R and the effective trap length / , of the correlation length l R in Eq. (459). The 
effective trap length / = 1 / ^/mwo and the effective trap frequency u are defined in Eq. (532). 
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More strictly, the system state, produced by the trap modulation, depends not merely on 
the modulation amplitude V mo d, but on the amount of the total energy pumped into the sys- 
tem, playing the role of the effective temperature in the nonequilibrium system [237]. For an 
alternating field, with the driving frequency u> and the related period t mo( i = 2tt/u, acting on 
the system during the total time of experiment t exp , the pumped energy is associated with the 
effective temperature (624). This energy should be treated as the effective modulation amplitude. 
Thus, instead of the localization length (460), we get 



2 



hoc — — n T/2 ;3 — ( T7^ ) ' (631) 



m 2 V 2 l 3 V V 

where, for clarity, the Planck constant is restored. 

Depending on the relation between the localization length (631) and the trap length Iq, there 
can exist the following states: 

hoc > lo (super fluid) , 

a < hoc < lo (heterophase fluid) , (632) 
hoc < a (chaotic fluid) . 

The superfluid state here includes all types of superfluids, the regular superfluid having no 
vortices, the vortex superfluid with a small number of vortices, and the turbulent superfluid 
with a random tangle of many vortices. This classification, in terms of the pumped energy, reads 
as follows: 

Vq < hu (superfluid) , 

hu <V < f\WQ y/lo/a (heterophase fluid) , (633) 
Vq > hujo \Jlo/a (chaotic fluid) . 

Chaotic fluid is a strongly fluctuating system having neither long-range order nor even local 
order. It resembles the state of weak turbulence [242] or the chaotic state [243]. Qualitatively, 
the overall scheme, representing the sequence of states arising under the action of an alternating 
field, with respect to the amount of the pumped energy Vq, is shown [237] in Figure 1. 

After the external modulation field is switched off, a finite quantum system relaxes to its 
equilibrium state during the relaxation time defined by particle collisions, the trap size, and trap 
shape [244]. The relaxation time becomes quite long for quasi-one-dimensional traps, where it 
may last, without equilibration for thousands of collisions between the oscillating Bose-condensed 
droplets [245] . This is because the one-dimensional system with local interactions is the integrable 
Lieb-Liniger system [246,247]. And quasi-one-dimensional systems, approaching integrability, 
display very long equilibration times. 

In a quasi-one-dimensional trap, collisions, restricted to the motion of particles in the axial 
direction, with the particles remaining in the same transverse ground state, are not accompanied 
by energy change, hence, do not lead to thermalization. For such two-body collisions, equili- 
bration and thermalization occurs only under transverse excitations. The corresponding rate of 
populating the radially excited modes by pairwise collisions, can be estimated from the Fermi 
golden rule that, at low temperature T < u±, gives [248,249] the rate 

T 2 « 2.8u ± (exp(-2u ± /T) , (634) 

where the dimensionless parameter 

C = Pid -f 
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is expressed through the three-dimensional scattering length a s , transverse oscillator length l±, 
and the one-dimensional density 

,2 N 
pid = pnl± = - . 

At temperature tending to zero, this rate is exponentially suppressed. However, it can be essential 
for finite temperatures. 

For quasi-one-dimensional traps at low temperatures, the three-body collision rate [249] can 
become important, 

T 3 « 6.9CV . (635) 

The states, described above, have been experimentally realized by means of the trap mod- 
ulation [241,250]. The whole diagram, showing the dependence of the produced states on the 
modulation amplitude and modulation time has been presented [250], starting from the regular 
superfluid, through vortex superfluid, to turbulent superfluid, and to heterophase fluid. 

12 Conclusions 

In this review, the basic theoretical problems have been considered, arising in the description 
of systems with Bose-Einstein condensate. The solutions to these problems are elucidated. The 
main conclusions can be briefly summarized as follows. 

(i) The global gauge symmetry breaking is the necessary and sufficient condition for Bose- 
Einstein condensation. This is an exact mathematical fact. The symmetry breaking results in 
the appearance of both, the condensate fraction and anomalous averages. The latter cannot be 
neglected without destroying the theory self-consistency. Omitting the anomalous averages is 
principally wrong, yielding unreliable and often unreasonable results. 

(ii) The Hohenberg-Martin dilemma of conserving versus gapless theories is resolved by intro- 
ducing two Lagrange multipliers guaranteeing the validity of two normalization conditions, for 
the numbers of condensed and uncondensed particles. The use of these two Lagrange multipliers 
is necessary as soon as the global gauge symmetry has been broken. 

(iii) Bose-Einstein condensed systems in strong spatially random potentials can be described 
by means of the method of stochastic decoupling. Perturbation theory with respect to the 
strength of disorder can fail, leading to incorrect conclusions. 

(iv) Thermodynamically anomalous fluctuations of any observable quantities are strictly pro- 
hibited in all equilibrium statistical systems, irrespectively of the used representative statistical 
ensemble. Thermodynamically anomalous particle fluctuations, of either condensed or uncon- 
densed particles, cannot exist in Bose-condensed systems. The occurrence of thermodynamically 
anomalous fluctuations can be due only to calculational mistakes. 

(v) The method has been suggested of generating nonground-state condensates of trapped 
particles. The method can be realized by applying alternating external fields modulating either 
the trapping potential or particle interactions. This makes it possible to create different types 
on nonground-state condensates, such as coherent modes, turbulent superfluids, and heterophase 
fluids. 
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Figure 1: Scheme of the sequence of states for a trapped Bose-condensed system subject to the 
action of an alternating external field, with the increasing pumped energy Vq. 
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